Geostatistical Evaluation of a Porphyry Copper Deposit Using Copulas

Abstract

Kriging has some problems such as ignoring sample values in giving weights to them, reducing dependence structure to a single covariance function, and facing negative confidence bounds. In view to these problems of kriging in this study to estimate Cu in the Iju porphyry Cu deposit in Iran, we used a convex linear combination of Archimedean copulas. To delineate the spatial dependence structure of Cu, the best Frank, Gumbel, and Clayton copula models were determined at different lags to fit with higher-order polynomials. The resulting Archimedean copulas were able to describe all kinds of spatial dependence structures, including asymmetric lower and upper tails. The copula and kriging methods were compared through a split-sample cross-validation test whereby the drill-hole data were divided into modeling and validation sets. The cross-validation showed better results for geostatistical estimation through copula than through kriging in terms of accuracy and precision. The mean of the validation set, which was 0.1218%, was estimated as 0.1278% and 0.1369% by the copula and kriging methods, respectively. The correlation coefficient between the estimated and measured values was higher for the copula method than for the kriging method. With 0.0143%² and 0.0162%² values, the mean square error was substantially smaller for copula than for kriging. A boxplot of the results demonstrated that the copula method was better in reproducing the Cu distribution and had fewer smoothing problems.

1. Introduction

Geostatistical estimation through kriging has become a standard method in mining engineering and earth sciences [1,2,3,4,5]. However, kriging weights given to samples are determined regarding the samples’ spatial configuration and a variable’s spatial continuity [6,7]. There are some problems with kriging, namely (1) it uses covariance, which describes spatial continuity by a single function imposing simplification in the estimation process [6,8], (2) it cannot give different weights to different sample values with the same spatial configurations [9], and (3) it provides symmetric confidence interval for the estimates, which may result in negative concentrations for small values [10]; however, the latter problem does not arise in nonlinear kriging.

The application of copulas in conjunction with geostatistics would be an excellent choice to solve the above-mentioned problems. A copula is a function that represents the joint distribution of variables that are uniformly distributed on [0, 1] [11]. Copulas have been used in different studies to evaluate a variable’s dependence structure through models such as student-t, Gaussian, chi-square, Gumbel, Clayton and Frank [12,13,14,15,16,17,18,19]. A spatial copula tackles the above-mentioned first and second problems by considering both sample amounts and spatial dependence structure to assign weights to the conditioning data. The confidence intervals generated by a spatial copula can take any shape and are not necessarily symmetric; therefore, appropriate estimates of the probability density function of the variable under study would result in reasonable confidence intervals.

Spatial copula was first introduced by Bárdossy [20] for geostatistical estimation of groundwater qualities based on the Gaussian and chi-square copulas. Because of its simple structure and ease of use, the Gaussian copula has been the topic of many studies [21,22]. However, like the other elliptical copulas such as the Student-t copula [23], the Gaussian copula cannot deal with variables with asymmetric spatial dependence structures. Therefore, applying elliptical copulas in modeling a variable’s spatial dependence structures could lead to unrealistic and simplified delineations [20,24,25]. The symmetric nature of Gaussian copula has forced researchers to use new copula models such as chi-square copulas [26]. However, these copulas are asymmetric with larger upper tails. Therefore, chi-square copulas and their variants cannot handle dependence structures with smaller lower tails [27]. To handle specific tails, some other studies proposed the Archimedean family of copulas such as the Gumbel, Clayton, Frank, and Joe models. However, none of these copulas can describe all possible structures. Therefore, the idea of combining them to get new Archimedean copulas that are not necessarily symmetric and theoretically can handle all types of asymmetric tails has become attractive [14,28,29,30]. Moreover, convex combinations of Archimedean copulas can be used with vine copula [8,31,32,33]. Copulas have provided flexible tools for describing a variable’s spatial dependence structure and for giving promising results in geostatistical estimations [10,34]. Therefore, in this study, due to the abovementioned advantages, the convex linear combination of Archimedean copulas (CLCAC) was used in the estimation of Cu in the Iju porphyry Cu mine, Iran. By running a split-sample cross-validation test, the results of geostatistical estimation through copulas was compared to results of kriging.

This paper is structured as follows. The theory of geostatistical estimation through copulas, the properties of Archimedean copula, and their convex linear combinations are presented in the Methods section. A case study including geological setting, data description, comparison of copula and kriging methods, and block modeling process are given in the Case Study section. The conclusions are presented in the last section.

2. Materials and Methods

2.1. Copula

A copula, $C$, is a function that can be used to delineate the dependence structure of variables, thus:

$$C:{\left[0,1\right]}^n\to \left[0,1\right]$$

(1)

If one of the n random variables takes zero, the copula function takes zero as well. According to Sklar [11], any n-variate distribution $F\left(Z_1, \dots ,Z_n\right)$ can be represented by its margins, $F_{Z_i}\left(Z_i\right)$, and n-dimensional $C$, thus:

$$F\left(Z_1, \dots , Z_n\right)=C\left(F_{Z_1}\left(Z_1\right),\dots ,F_{Z_n}\left(Z_n\right)\right)$$

(2)

The transformation of a variable to standard uniform distribution makes the copula independent of the margins. Standard uniform margins can be achieved through a simple histogram transformation. For continuous margins, a copula is unique and can be written in terms of mutual dependence structure regardless of margins. Therefore, the density of a copula, $c$, can be achieved from the following derivative:

$$c\left(u_1, \dots , u_n\right)=\frac{{\partial }^{n}C\left(u_1,\dots ,u_n\right)}{\partial u_1\dots \partial u_n},$$

(3)

and the conditional copula can be calculated as:

$$C\left(u_1 U_2=u_2, \dots , U_n=u_n\right)=\frac{{\partial }^{n-1}C\left(u_1,\dots , u_n\right)}{\partial u_2\dots \partial u_n}\times \frac{1}{c\left(u_2, \dots ,u_n\right)},$$

(4)

2.2. Spatial Copula

Assume that $Z$ is a second-order stationary random field sampled at $Ɲ$ points. Second-order stationarity provides the following property for two sets of samples separated by a vector, $\mathit{h}$:

$$P\left(Z\left(x_1\right)<v_1, \dots , Z\left(x_k\right)<v_k\right)=P\left(Z\left(x_1+\mathit{h}\right)<v_1, \dots , Z\left(x_k+\mathit{h}\right)<v_k\right)$$

(5)

where $v_1,\dots ,v_k$ are some possible values and $Z\left(x_1\right)$ is the value of $Z$ at location $x_1$. For a single variable, the bivariate representation of a spatial copula at vector, $\mathit{h}$, takes the following form:

$$C_s\left(\mathit{h},u,u\right)=P\left(F_Z\left(Z\left(x\right)\right)<u, F_Z\left(Z\left(x+\mathit{h}\right)\right)<u\right)=C\left(F_Z\left(Z\left(x\right)\right), F_Z\left(Z\left(x+\mathit{h}\right)\right)\right)$$

(6)

2.3. Copula Modeling

The empirical copula of a variable can be achieved as:

$$\mathit{Y}\left(\mathit{h}\right)=\left\{F_n\left(Z\left(x_i\right),Z\left(x_j\right)\right)| \left(x_i-x_j\right)\approx \mathit{h} or \left(x_j-x_i\right)\approx \mathit{h}\right\}$$

(7)

where $F_n\left(z\right)$ is the empirical distribution function obtained from observations $z_1,\dots ,z_{Ɲ}$. Various theoretical copulas,$C^{*}$, can be fitted to the empirical copula of Equation (7) to find the best function for which the null hypothesis at a significance level, α, is not rejected. The procedure can be implemented through a two-sample Kolmogorov-Smirnov test for equality of functions [20], thus:

$$D_{K-S}=\sup \left\{\left|C^{*}\left(u_1,u_2\right)-C\left(u_1,u_2\right)\right| \left(u_1,u_2\right)\in {\left[0,1\right]}^2\right\}$$

(8)

where $u_1$ and $u_2$ are the measured values of the two samples after standard uniform transformation.

Fitting the copula associated with different lag separation vectors is equivalent to fitting the bivariate distributions of the random field for these separation vectors. In general, the fit does not guarantee the existence of a random field associated with such bivariate distributions so that there may be a problem of internal consistency of the copula model, which has been pointed out by Matheron [35] and has still not received an answer. This is a counterpart of attempting to avoid a simplified spatial continuity model.

2.4. Combination of Archimedean Copulas

Assume that $C\left(u_1,u_2\right)$ is an Archimedean copula presented as:

$$C\left(u_1,u_2\right)={\phi }^{-1}\left(\phi \left(u_1\right)+\phi \left(u_2\right)\right),0\le u_1, u_2\le 1$$

(9)

For Archimedean copulas, $\phi$ is the generator that is a convex decreasing function on $\left[0,1\right]$ with $\phi \left(0\right)=\infty$ and $\phi \left(1\right)=0$ [36,37,38]. The inverse function of $\phi$ is ${\phi }^{-1}:\left[0,\infty \right]\to \left[0,1\right]$ for which ${\phi }^{-1}\left(0\right)=1$ and ${\phi }^{-1}(\infty )=0$. Any function with these mentioned properties can be used as a generator.

Table 1. Common Archimedean copulas and their properties.

Copula Name	Generator ${\mathit{\phi }}_{\mathit{\theta }}\left(\mathit{t}\right)$	Parameter Range	$\mathbf{Copula} {\mathit{C}}_{\mathit{\theta }}\left({\mathit{u}}_1,{\mathit{u}}_2\right)$	${\mathit{\lambda }}_{\mathit{L}}$	${\mathit{\lambda }}_{\mathit{U}}$
Gumbel $C^G\left(u_1,u_2\right)$	${(-lnt)}^{\theta }$	$\theta \ge 1$	$exp\left\{-{\left[{(-ln{u}_1)}^{\theta }+{(-ln{u}_2)}^{\theta }\right]}^{\frac{1}{\theta }}\right\}$	0	$2-2^{\left(\frac{1}{\theta }\right)}$
Clayton $C^C\left(u_1,u_2\right)$	$t^{-\theta }-1$	$\theta \in \left[-1,\infty \right)\backslash \left\{0\right\}$	${\left[u_1{}^{-\theta }+u_2{}^{-\theta }-1\right]}^{-\frac{1}{Ɲ}}$	$2^{-\frac{1}{\theta }}$	0
Frank $C^F\left(u_1,u_2\right)$	$-ln\left(\frac{exp\left(-\theta t\right)-1}{exp\left(-\theta \right)-1}\right)$	$\theta \in \mathbb{R}\backslash \left\{0\right\}$	$-\frac{1}{\theta }ln\left[1+\frac{\left(exp\left(-\theta u_1\right)-1\right)\left(exp\left(-\theta u_2\right)-1\right)}{exp\left(-\theta \right)-1}\right]$	0	0

presents some common Archimedean copulas such as the Gumbel, Clayton, and Frank copulas and their properties. Each copula generator has one parameter, $\theta$, defined in a specific range. The lower and upper tails of copulas are denoted by ${\lambda }_L$ and ${\lambda }_U$, respectively. For a bivariate distribution, lower and upper tail dependences refer to the degree of dependence in the corner of the lower-left quadrant and upper-right quadrant, respectively.

It can be seen from Table 1 that the Gumbel and Clayton copulas are asymmetric with upper and lower tails, respectively. None of these copulas is capable of describing all kinds of spatial dependence structures on its own. However, the CLCAC provides a flexible tool that can be fitted to an extensive range of structures [39,40]. The $C^{GC}$, given as the following formula, combines the copulas, thus:

$$C^{GC}=\{\begin{array}{c}{\alpha }_1{C}_{h_1}^G\left(u_1,u_2\right)+\left(1-{\alpha }_1\right)C_{h_1}^C\left(u_1,u_2\right)\hspace{1em}\hspace{1em}if\hspace{1em}\hspace{1em}h=h_1\\ \begin{array}{c}{\alpha }_2{C}_{h_2}^G\left(u_1,u_2\right)+\left(1-{\alpha }_2\right)C_{h_2}^C\left(u_1,u_2\right)\hspace{1em}\hspace{1em}if\hspace{1em}\hspace{1em}h=h_2\\ ⋮\\ {\alpha }_l{C}_{h_l}^G\left(u_1,u_2\right)+\left(1-{\alpha }_l\right)C_{h_l}^C\left(u_1,u_2\right)\hspace{1em}\hspace{1em}if\hspace{1em}\hspace{1em}h=h_l\end{array}\end{array}$$

(10)

where $C_{h_1}^G$ and $C_{h_1}^C$ stand for the best Gumbel and Clayton copulas, respectively, that are fitted to the empirical copulas at lag distance $h_1$, and ${\alpha }_i (i=1,2,\dots ,l$) is the weight given to the first copula in the combination. Weights assigned to each copula in the combinations can be found by testing a series of values in $\left[0,1\right]$ whereby a specific step size is applied or they can be found by defining an objective function based on Kolmogorov-Smirnov test and optimizing it using methods such as simulated annealing. By calculating the derivative of Equation (10), the copula density and conditional density at the estimation points can be obtained. By integrating the conditional density and calculating its median value, 2.5% lower bound, and 97.5% upper bound, the estimated value and 95% confidence interval can be obtained.

2.5. Steps of Prediction

After transforming data into standard uniform distribution through $\mathrm{U}\left({\mathrm{x}}_{\mathrm{i}}\right)={\mathrm{F}}_{\mathrm{Z}}\left(\mathrm{Z}\left({\mathrm{x}}_{\mathrm{i}}\right)\right)={\mathrm{u}}_{\mathrm{i}}$ and calculating CLCAC, the conditional copula density, ${\mathrm{c}}_{{\alpha }_i}^{\mathrm{GC}}\left({\mathrm{u}}_0{|\mathrm{u}}_{\mathrm{i}}\right)$, can be achieved at the prediction point, ${\mathrm{x}}_0$. Weights given to the copulas, ${\alpha }_i$, depend on the separation distance between the prediction location, ${\mathrm{x}}_0$, and sample locations, ${\mathrm{x}}_{\mathrm{i}}$. Prediction through CLCAC can be performed as follows:

(i)

Transform the variable into uniform distribution through $U=F_Z\left(Z\right)$.

(ii)

Choose lag distances and lag tolerance for calculating empirical copulas for them.

(iii)

For each lag, estimate the best copula parameters and find their best combination.

(iv)

To calculate copula density at distances other than lag distances, model copula parameters and weights given to them.

(v)

For each conditioning sample, put its distance from the prediction point into the fitted models of the previous step to get its density at the prediction point.

(vi)

Compute conditional cumulative distribution function at the prediction point in order to get weights given to the samples. Multiply sample values by their weights, sum up the multiplication results and back-transform them into original data space.

3. Case Study

3.1. Geology

The Iju area, within 30°31′45″ N to 30°33′05″ N and 54°56′10″ E to 54°57′30″ E (Figure 1), is located 42 km NW of Shahre-Babak county, Kerman province, and 140 km NW of the Sarcheshme Cu mine. The Iju deposit is situated on the SE part of the Urmia–Dokhtar magmatic belt [41,42,43], which is characterized by numerous Cu deposits such as the Chah-Firouze, Sarcheshme, Meiduk, and Chah-Messi (Figure 1).

The Iju deposit is situated in a mountainous area with Eocene–Paleocene pyroclastic volcanic rocks intruded by Miocene quartz-diorite and tonalite [44] (Figure 1). Based on field investigations and core sample analyses, quartz diorite and tonalite are respectively intruded into the host rocks (Figure 1). Then, alteration and Cu mineralization has occurred. Extensive phyllic alteration and surface occurrences of malachite, chalcanthite, jarosite, and iron hydroxides minerals are evidence of Cu mineralization in the study area (particularly in the southern part). Phyllic alteration with sericite, quartz, pyrite, and chlorite minerals is widespread and detected in almost all of the drill holes in the central part of the study area. In contrast, potassic alteration, characterized by the formation of secondary-biotite, and K-feldspar veins, is limited and only observed in some drill holes in the southern part. Propylitic alteration is widespread in the host volcanic and pyroclastic rocks. Argillic alteration is detectable on and near the surface.

Figure 1. Location and geological maps of the Iju deposit [45].

Figure 1. Location and geological maps of the Iju deposit [45].

In the Iju deposit, Cu mineralization exists in the form of disseminations and stockworks. Chalcopyrite is the main Cu mineral, often seen along with pyrite. Magnetite is observed as veins and veinlets in the central part of the study area. Tiny amounts of gypsum, anhydrite and molybdenite also exist in the study area. Due to the potassic, propylitic, and extensive phyllic alterations, and the mineralization style, the Iju is classified as a porphyry Cu deposit [44,46]. Four hydrothermal mineralization zones, namely leached, oxidized, supergene, and hypogene, are found in the area. Thicknesses of the leached and supergene zones vary between 10 and 50 m and between 5 and 50 m, respectively. There are two high-grade Cu zones in the northern and southern parts of the study area, joining with a low-grade central part. As shown in Figure 2, the main portion of the Iju resource is associated with the hypogene zone, and this study was focused on this zone.

3.2. Data Description

Core samples from 39 drill holes (Figure 3) were used in the estimation process. The average distance between adjacent drill holes and the average core length were 50 m and 2 m, respectively. Therefore, 5458 2 m composites with histogram and summary statistics, respectively shown in Figure 4 and

Table 2. Summary statistics for 2-m composite values.

Variable	Number	Mean	Variance	Skewness	Minimum	Maximum
Cu	5458	0.1680	0.0250	2.37	0.0006	1.3884

, were generated for further analysis. Due to strong deviation from normal distribution, the data were first transformed into normal scores in the SGeMS program. Then, the spatial dependence structure of the normal scores of Cu was assessed in the same program by plotting experimental variograms in the downhole and different horizontal directions. The dependence structure of Cu did not show any zonal anisotropy with variograms reaching the sill value of 1. The Cu data showed geometric anisotropy with a maximum variogram range of 260 m in the downhole direction. Directional horizontal variograms demonstrated almost identical ranges (175 m). Variograms of Cu were fitted with a nested model, including a nugget effect, a short-range exponential structure, and a long-range spherical structure. Figure 5 and

Table 3. Model variogram parameters.

Variogram	Nugget	Structure #1	Structure #2	C1	C2	Range 1 (m)	Range 2 (m)
Horizontal	0.1	Exp	Sph	0.15	0.75	7	175
Downhole						10	260

, respectively, present the experimental variograms and fitted models, and model variogram parameters used in the kriging method.

Spatial dependence evaluation through copulas was performed using a MATLAB code developed by the first author. An omnidirectional evaluation was performed at 30 lags considering lag distances and tolerance of 2 m and 1 m, respectively. Lagged scatterplots for the variable were obtained at different distances and some of them are shown in Figure 6. For each lag distance, the best Archimedean copula parameters $(\theta )$, which appropriately delineate the dependence structure, were obtained (Figure 7). Moreover, various $\alpha$ weights, between 0 and 1, considering a specific step size, were given to each copula in the Gumbel–Clayton, Clayton–Frank, and Gumbel–Frank combinations to find their best portions (Figure 8). The empirical points in Figure 8 were fitted by appropriate polynomial functions to calculate copula parameters at other lags. From Figure 8, it can be seen that, at 60 m, the $\alpha$ weight given to the Clayton copula in the Clayton–Gumbel combination starts from 0.1 and reaches 0.65. Therefore, the spatial dependence structure of Cu had a strong positive tail at close distances and relatively strong negative tails at large lags. In some distances, especially when h increases to infinity, the fitted model may result in unreasonable $\alpha$ values. For such cases, the developed MATLAB code replaces negative values by zero and changes large amounts to 1.

3.3. Split-Sample Cross-Validation

Comparison of kriging and copula methods was performed by running a cross-validation test based by dividing the data into two sets: (1) a modeling set with 4237 samples and (2) a validation set with 1221 points. Variograms and copula dependence structures were evaluated based on the modeling set. Then, the obtained parameters were used in the estimation process at the validation points. Estimation through both methods was performed using the same number of conditioning data (17 samples) to obtain comparable results. Ordinary kriging was performed twice using the original values and also based on normal scores. As expected, the normal score-based kriging performed substantially better than kriging based on the original values. Therefore, hereafter, we only demonstrate the results of the normal score-based kriging. The estimated values were compared to the measured values to assess the accuracy and precision of the estimators (

Table 4. Cross-validation results.

Estimation Method	Measured Mean	Estimated Mean	Mean Squared Error (Perfect Value Is 0)	Correlation (Perfect Value Is 1)
Kriging	0.1218	0.1369	0.0162	0.42
Copula		0.1278	0.0143	0.51

). The copula method outperformed kriging based on having smaller mean squared error, better mean reproduction, and a higher correlation coefficient between the estimated and measured values. Because of having closer quartile values, the data distribution was reproduced better in the copula estimates (Figure 9). The cross-validation results indicate the importance of estimation through copulas for a highly skewed data such as the Cu data. Kriging was significantly affected by a few large sample values and the L-shape distribution, and it showed a tendency for over-estimation.

3.4. Estimation and Results

After the cross-validation test, the best copula and variogram parameters were used in the block estimation process. A 3D geological model of the hypogene zone and a block model including 16,724 blocks ($25 \mathrm{m}\times 25 \mathrm{m}\times 12.5 \mathrm{m}$) were created (Figure 2). To reduce the change of support problem, estimation was performed by dividing the main blocks into 27 equal size sub-blocks, running estimation at center points, and taking the average of the estimated values [47,48,49]. The estimated block models, box plots of the block estimates, and summary statistics of the results are, respectively, given in Figure 10 and Figure 11, and

Table 5. Summary statistics for the data and the kriging and copula estimates.

	Mean	First Quartile	Median	Third Quartile	Variance	Min	Max
Data	0.168	0.060	0.125	0.222	0.025	0.0006	1.388
Kriging	0.134	0.052	0.102	0.191	0.011	0.008	0.913
Copula	0.144	0.047	0.104	0.200	0.019	0.006	1.099

.

Except for the first quartile, other statistical properties of data such as the mean, median, third quartile, minimum, maximum, and variance values were reproduced better in the copula estimates than in the kriging estimates. The variance of the kriging estimates was 0.011 (%²), which was less than half of the data variance of 0.025. From this aspect, the copula results with 0.019 (%²) variance showed potential for reducing the smoothing problem of kriging. This issue is reflected by the minimum and maximum values of the copula estimates, which were smaller and larger compared to the minimum and maximum values of the kriged estimates, respectively.

Box plots of the estimation errors calculated for the results of both methods demonstrated significantly smaller errors for the copula method (Figure 12). The mean of estimation errors for the copula method was 1.80 compared to 4.82 of the estimation errors for the kriging method. Except for a limited number of outliers, the estimation errors of the copula method were lower than 3.13; however, 82% of the kriging estimation errors were larger than 3.38.

4. Conclusions

The problems of kriging (e.g., simplifying the dependence structure, lacking the ability to consider sample values in assigning weights to them, and the possibility of producing absurd confidence intervals) motivate the use of new methods such as spatial copula in geostatistical applications. Therefore, in this study, the CLCAC was used to estimate the copper grade in the Iju porphyry copper deposit and it was compared with kriging. Due to the flexible tails of the CLCAC, the proposed method delineated successfully the spatial dependence structure of the variable under study (namely Cu data). Comparison of the kriging and copula estimates, through a split-sample cross-validation test, showed advantages of the proposed method over kriging. The copula estimates had lower mean squared error, higher correlation coefficient with data, better reproduction of mean, and closer distribution to that data. The same advantages were also be seen in the block estimation results of the copula and kriging methods. The copula estimation errors were significantly smaller than the kriging estimation errors. Moreover, the copula estimates had higher variance than the kriging estimates; thus, using a spatial copula mitigates the smoothing problem of kriging.

The results of this study suggest that the copula approach is dramatically outperforming ordinary kriging. This issue can be a consequence of choosing the right dataset with strong deviation from normal distribution and testing an advanced method against a method (kriging) that is too simple to deal with the special challenges of that dataset.

In this study, univariate estimation was performed based on two-point statistics. However, it is known that methods relying on two-point information can be outperformed by high order methods such as high-dimensional copulas, which use high order dependences. Therefore, the application of multivariate cases and high-dimensional copulas remains as a future study.

Another point that remains as a future study is associated with integrating directional anisotropy with the method. The copula approach applied in this study uses an omnidirectional model. However, one can upgrade the method by dividing a 3D space into some equal size slices. Therefore, the spatial continuity of the variable under study can be described directionally in each slice.