Integration of Dual Border Effects in Resource Estimation: A Cokriging Practice on a Copper Porphyry Deposit

: Hierarchical or cascade resource estimation is a very common practice when building a geological block model in metalliferous deposits. One option for this is to model the geological domains by indicator kriging and then to estimate (by kriging) the grade of interest within the built geodomains. There are three problems regarding this. The ﬁrst is that sometimes the molded geological domains are spotty and fragmented and, thus, far from the geological interpretation. The second is that the resulting estimated grades highly suffer from a smoothing effect. The third is related to the border effect of the continuous variable across the boundary of geological domains. The latter means that the ﬁnal block model of the grade shows a very abrupt transition when crossing the border of two adjacent geological domains. This characteristic of the border effect may not be always true, and it is plausible that some of the variables show smooth or soft boundaries. The case is even more complicated when there is a mixture of hard and soft boundaries. A solution is provided in this paper to employ a cokriging paradigm for jointly modeling grade and geological domains. The results of modeling the copper in an Iranian copper porphyry deposit through the proposed approach illustrates that the method is not only capable of handling the mixture of hard and soft boundaries, but it also produces models that are less inﬂuenced by the smoothing effect. These results are compared to an independent kriging, where each variable is modeled separately, irrespective of the inﬂuence of geological domains. produces higher recovery functions, particularly higher mean grade, for further economic evaluations. For instance, the fraction of tonnage, mean grade, and metal quantity above 0.3% Cu for Case I and Case II are: 8.4% and 3.4%, 0.453% and 0.407%, and 0.038% and 0.014%, respectively.


Introduction
Mineral resource estimation is the first stage in the exploration and development of a mining project [1][2][3], and its reliability in the downstream of mining activities such as long-and short-term mine planning is inevitable [4,5]. The traditional approach in metalliferous deposits is based on first identifying the layout of the geological domain in the region and then estimating the grade of interest within each geodomain, which is called a hierarchical or cascade approach [6][7][8][9]. The former can be implemented by any implicit or explicit techniques, and the latter by using any interpolation or geostatistical techniques. To name a few, radial basis function (RBF) [10,11] and wireframing [12][13][14] are two common approaches for implicit and explicit geological modeling, respectively. An alternative to explicit modeling is indicator kriging [15,16], where an indicator can characterize each domain; afterwards, through a conditional cumulative probability function, the most probable domain can be estimated at unsampled locations. Pereira et al. (2017) [17] showed that indicator kriging produces better results compared to wireframing, which is very time-consuming. However, this simple-to-use method sometimes produces very patchy and unstructured maps, which makes it difficult to interpret in further steps of an ore body evaluation. An important characteristic in modeling of grade, given the geological information, is the condition of grade variability throughout the geological boundaries. A boundary is soft if the grade shows gradual changes once crossing the adjacent geological domains, and is hard if the changes across the geological boundary are abrupt. Another issue in this hierarchical technique of resource estimation is yielding hard boundaries in the continuous variable when crossing the border of geological domains. This means that this cascade workflow neglects the importance of the border effect if there is a soft boundary in the contact of two adjacent domains [18,19]. To come up with this problem, Ortiz and Emery (2006) [20] presented an algorithm based on the incorporation of detailed information about the adjacent domains up to a given radius by using ordinary kriging when estimating the main variable in a copper deposit. The method showed satisfying results in the case of having pure soft boundaries when compared to other resource estimation techniques; however, it is not able to handle the hard boundaries. In practice, complications arise when one faces a mixture of hard and soft boundaries in a deposit. An alternative is to use indicator kriging to model first the probability of estimation domains in the region and use a separate kriging to model the grade of interest. The final grade estimate can be achieved through weighting the antecedent estimates in accordance with the probability of pertaining to the domain of interest [21,22]. This approach is suitable for hard boundaries, but it is also effective, with less precision, for the soft boundary. This method is presented for the case where domains are identified through partitioning the grade (i.e., grade domains), but can still be used for situations where estimation domains are characterized by geological information. One significant inefficiency in this method is that, since indicator kriging is used independently to model the estimation domains, there is a similar issue of patchy and unstructured boundaries of some geodomains in the resulting maps. Kasmaee et al. (2019) [23] presented the alternative of a transitional boundary, where there is an overlapping of geological domains. Although the method diminishes the magnitude of the disconnections in the domain boundaries, it does not eliminate them entirely [24].

Resource Modeling
In mineral resource modeling, one is usually dealing with the grade of interest in a deposit, which is a continuous variable, and geological information, which is a categorical variable. The continuous variables are numeric variables that can have an infinite number of values between any two values, while the categorical variables are the variables that can take one of limited natural numbers, assigning each individual to particular geological information. The categorical variable can be converted to indicator variables by Equation (1).
Hierarchical/cascade resource modeling is a typical practice in the process of building a geological block model. In this method, first the geological attributes are modeled by any implicit or explicit techniques, and the continuous variable(s) are then estimated inside the already-built envelope of the geological domains, taking into account the data that only belong to that corresponding domain [2,3]. The definition of these domains can be stemmed based on the geological variability, such as alterations, rock types, and mineralized zones, in most of the metalliferous deposits. In nonmetalliferous deposits such as coal, however, this might be different. In these deposits, seam layers are the main estimation domains. Another approach for deriving these categorical domains is truncating the main continuous variables into several domains based on the prespecified cut-offs. In this technique, the variable of interest is separated into several mineral classes such as low-, medium-, and high-grade domains. However, the resulting domains in this method should be compliant with geological interpretations of the deposit.
The hierarchical method is widely used in commercial software programs. One restriction for this approach is that the behavior of continuous variable(s) across the geological boundaries should be hard (i.e., abrupt transition), although this characteristic is overlooked in most of the modeling workflow. One reason for such independent modeling of the continuous variable(s) inside each domain is related to the geological character-istics of the main variable(s), which make it distinguishable in each different geological domain. In other words, separate modeling of continuous variable inside each domain guarantees the homogeneity of the modeled continuous variable inside each domain. This might be true if the behavior of the continuous variable across the geological boundary is hard, but it is not true if the transition of the continuous variable across the boundary is soft, and there is a mixture of hard and soft boundaries. Regarding this latter case, an algorithm was used in this study based on the cokriging of continuous and categorical/indicator variables simultaneously, which can produce a mixture of border effects in the final estimation results.
The algorithm for one continuous variable and n indicator variables is as follows: • Converting the geological domain (i.e., categorical data) into indicators at sample locations x: In the case of having n categorical data, the sum of the n indicator variables is one at sample locations.

•
Computation of experimental direct and cross-variograms or direct and cross-covariances among the continuous variables and indicators. A linear model of coregionalization can be used for fitting purposes [25]. In this step, modeling a coregionalization calls for inferring (n + 1) × (n + 2)/2 direct and cross-variograms or direct and cross-covariances. Assuming numerous indicator variables, a semiautomatic iterative algorithm [26] can be used for fitting purposes. • Using a cokriging system to estimate the n + 1 variables. The "1" here refers to one continuous variable. If the mean of all n + 1 variables is known, then a simple cokriging is used, and if the mean is unknown, an ordinary cokriging can be considered. The neighborhood can also be moving or unique, depending on the size of conditioning data. The implementation of this cokriging system is provided in the following section.
Since the results of indicator cokriging suffer from an order relation problem [26] for the resulting estimated indicators, a postprocessing step is needed to correct this issue.

The Cokriging System
The cokriging system in this algorithm is dealing with one continuous variable and n indicators, considering their joint spatial correlation structure. Provided that these variables are represented by second-order stationary random fields and that the n + 1 variable means are known, the simple cokriging estimate and estimation variance for the continuous variable (hereafter denoted with "c"), given n indicator variables (hereafter denoted with "i"), can be defined as below: where λ c α and λ i α are the weights allocated to the data Z c (x c,α ) and Z i (x i,α ), respectively, at the α − th data location, conveying either continuous or indicator variables (or both), x 0 is the location targeted for prediction, m c and m i are the mean values of continuous and indicator variables, and C c and C i/c are the direct and cross-covariances, respectively, between continuous and indicator variables.
The weights λ c α and λ i α required in Equations (2) and (3) can be achieved by means of solving the following system of equations: This system of cokriging system is based on one continuous and several indicator data. The indicator part of Equations (2)-(4) implies the traditional indicator cokriging system that now is integrated with including one continuous data in the estimation process. Hence, the resulting estimated values for indicators Z * i at target locations x 0 suffer from an order relation problem. Due to some negative weights of the cokriging system, one might face some negative values or very high values for indicators, which are more than 1 [15,21]. In this case, one needs to take into account two postprocessing steps in order to correct the estimated values of the indicators. The first action is to use bounds correction, meaning that all negative estimated values are set to 0 and all high values (more than 1) are set to 1. The second difficulty is that the summation of all estimated values may not reach 1. Regarding this, a normalization of obtained probabilities is needed [27]. The negative estimates also can be seen in the estimated values of continuous variable. One solution for this is to set these values to the minimum value of the original data.
Once the order relation problem is resolved, n + 1 estimated values are ready at target grid nodes. One of them is a continuous variable that can be used for resource estimation activities. In order to obtain the categorical maps, one only needs to calculate the most probable estimated values at target location and assign the corresponding category, exactly the same as the procedure applicable for finalizing the indicator cokriging. In the end, each target location reports two values-one is a continuous variable, and another is a deterministic geological domain. The proposed algorithm is shown through a case study.

Geological Setting
The Urumieh-Dokhtar magmatic belt (UMDB) extends in the NW to SE direction for over 1700 km and has a 4-50 km width. The current distance from the Zagrus suture is 50-80 km [28,29] (Figure 1). The magmatic belt is built over a basement comprising (I) Eocene-Oligocene magmatic rocks with a geochemical signature typical of a continental arc under extension [30] from~55 Ma to~37 Ma [31], (II) the Oligocene OIB-like magmatic rocks, indicating a transition from arc-like to "OIB-like", (III) the late Oligocene-Miocene arc magmatism [32], which indicates the beginning of collision, and (IV) the late Miocene-Quaternary magmatism, which corresponded with the onset of closure of the Tethys [33,34]. However, it should be mentioned that middle Eocene magmatism is volumetrically dominant over the others [31]. From a geographical perspective, the UDMB is divided into three segments. One is the northwestern part, including the Sungun, Haft Cheshmeh, Sunajeel, Niaz, Miveh Rud, Kighal, Ali Javad, and Saheb Divan domains, typically associated with the Oligocene-Miocene fertile intrusions. The second is the central part, situated among Saveh and Yazd, which contains the Qom-Kashan-Natanz, Nain, and Nain-Yazd areas ( Figure 1). This part comprises late Neoproterozoic to Paleozoic sedimentary rocks overlain by Mesozoic rocks. All these units were covered by volcanic rocks during the Paleogene and imposed on by mafic and felsic plutons from Oligocene to Miocene times [35]. The third is the southeastern part, famous as the Dehaj-Sarduiyeh magmatic belt (DSMB) (Figure 2), which covers the Anar-Sirjan-Kerman, Bam-Jiroft, and Bazman areas. The DSMB is characterized by huge calc-alkaline magmatic rocks containing late Eocene-Oligocene and) middle to late Miocene intrusions with barren and fertile affinities, respectively [36]. mafic and felsic plutons from Oligocene to Miocene times [35]. The third is the southeastern part, famous as the Dehaj-Sarduiyeh magmatic belt (DSMB) (Figure 2), which covers the Anar-Sirjan-Kerman, Bam-Jiroft, and Bazman areas. The DSMB is characterized by huge calc-alkaline magmatic rocks containing late Eocene-Oligocene and) middle to late Miocene intrusions with barren and fertile affinities, respectively [36].  The felsic intrusions are subdivided into granodiorites, granites, diorites, and quartz diorite porphyries [38] ( Figure 3A). Granodiorite is chiefly intruded into the Eocene volcanic rocks ( Figure 4) and shows a microgranular porphyritic texture with a semi-blocky, brittle morphology. Quartz (up to 40 vol%) and plagioclase (~30 vol%), and amphibole and biotite (~15 vol%) constitute main minerals whereas apatite, titanite, and iron oxide minerals are accessory phases. Plagioclase phenocrysts (~3.5 mm) are slightly to moderately altered to sericite, calcite, and hematite. The hornblende and biotite have been also replaced by calcite, chlorite, epidote, and iron oxide minerals (Figure 5c,d). The Bondar Hanza copper mineralization at the southeast of the DSMB is one of the southeastward [28] continuations of the UDMB. The main lithologies in the study area include the lower Eocene intermediate volcanic and volcano-sedimentary rocks, which are intruded by late Eocene-early Oligocene felsic intrusions ( Figure 2).
Quartz (up to 40 vol%) and plagioclase (~30 vol%), and amphibole and biotite (~15 vol%) constitute main minerals whereas apatite, titanite, and iron oxide minerals are accessory phases. Plagioclase phenocrysts (~3.5 mm) are slightly to moderately altered to sericite, calcite, and hematite. The hornblende and biotite have been also replaced by calcite, chlorite, epidote, and iron oxide minerals (Figure 5c,d).    The diorite and quartz diorite and other shared intrusive rocks are the main hosts, covering an area of over 1 km 2 , and show a micro granophyric texture ( Figure 3A). Petrographically, diorite contains laths of plagioclase (>90 vol%) and insignificant amounts of quartz, K-feldspar, and iron oxide minerals. Plagioclase (2-3 mm) has been altered to sericite and argillite minerals. The andesite lava, tuff, lithic tuff, and ignimbrite make up the main rock types of the Bondar Hanza volcanic sequence. Andesite is located in the southern part of Bondar Hanza and is commonly characterized by porphyritic texture with feldspars, biotite, and amphibole phenocrysts set in fine-grained groundmass. The subhedral plagioclase phenocrysts (>30 vol%) look altered, mainly along microfractures.  The diorite and quartz diorite and other shared intrusive rocks are the main hosts, covering an area of over 1 km 2 , and show a micro granophyric texture ( Figure 3A). Petrographically, diorite contains laths of plagioclase (>90 vol%) and insignificant amounts of quartz, K-feldspar, and iron oxide minerals. Plagioclase (2-3 mm) has been altered to sericite and argillite minerals. The andesite lava, tuff, lithic tuff, and ignimbrite make up the main rock types of the Bondar Hanza volcanic sequence. Andesite is located in the southern part of Bondar Hanza and is commonly characterized by porphyritic texture with feldspars, biotite, and amphibole pheno-

Presentation of the Dataset
The dataset belongs to the Bondar Hanza Copper porphyry deposit diamond drill hole campaign. Through the sample locations, one continuous variable, copper, and three categorical variables, namely rock type, alternation, and zone, are identified. In order to preserve the confidentiality of the dataset, the copper is multiplied by a scale factor. The former continuous variable is measured through chemical assaying, and the latter categorical variables are recorded by geological logging of the core samples. Initially, five alteration types (argillic (ARG); phyllic (PHY); propylitic (PRP); silicification alteration (SLC); potassic (POT)), seven rock types (andesite (AND); ignimbrite (IBG); lithic tuff (LTT); tuff (TUF); quartz diorite (QDI); granodiorite (GRD); metavolcanic rocks (MTV)), and seven different zones (hypogene (HYP); leaching (LEA); oxide (OXI); oxide-hypogene (OXI-HYP); oxide-supergene (OXI-SUP); supergene (SUP); supergene-hypogene (SUP-HYP)) were recognized. However, in order to start the process of modeling, we decided to group these variables so as to better manage them into the 3D resource modeling workflow. To do so, the number of samples in a category and their location in the deposit were two important components for dimension reduction. Table 1 shows how the grouping was implemented over the sample dataset. As can be seen, the dimensions of alteration, rock type, and zone were reduced from five to three, from seven to three, and from seven to two, respectively.
The location maps of the variables are shown in Figure 6. As can be observed, the copper is mostly distributed in the center of the deposit. The rock group C (RG-C), the alteration group C (AG-C), and the zone group A (ZG-A) are the most dominant rock group, alteration, and zone in the region, respectively. This was also verified by calculation of the proportion of each category in Table 2.
Minerals 2021, 11, x FOR PEER REVIEW 10 of 23 lar, no declustering was applied on this dataset. One assumes that the provided statistical parameters are representative and are not influenced by excessive and clustered data in any area. To model the copper in this deposit, one needs to identify the estimation domains, since the continuous variable is deemed homogeneous in each domain, where it can be modeled separately [3]. For this purpose, we needed to build the boundary of the geological domain first and then estimate the copper within the built geological domain. However, in this deposit, one continuous and three categorical variables are involved; therefore, each category group can be considered for the identification of the estimation domain. It is commonly proposed that the barren to sub-fertile intrusions are diorite and granodiorite [36,40,41], while fertile intrusions are mostly quartz monzonite, tonalite, and minor granodiorite [42][43][44]. Both granodiorite and granite are considered as the hosts of copper mineralization in the Bondar Hanza ( Figure 3A), in agreement with sub-fertile mineralization.
Four alteration stages are recognized in the Bondar Hanza copper deposit based on the paragenesis and texture relationships, including (i) potassic, (ii) phyllic, (iii) propylitic, and (iv) argillic alteration zones (Figure 3). Potassic alteration is considered by the existence of biotite, quartz, The statistical parameters for the copper were calculated, and the results are presented in Table 2. Since the sampling pattern is almost regular, no declustering was applied on this dataset. One assumes that the provided statistical parameters are representative and are not influenced by excessive and clustered data in any area. To model the copper in this deposit, one needs to identify the estimation domains, since the continuous variable is deemed homogeneous in each domain, where it can be modeled separately [3]. For this purpose, we needed to build the boundary of the geological domain first and then estimate the copper within the built geological domain. However, in this deposit, one continuous and three categorical variables are involved; therefore, each category group can be considered for the identification of the estimation domain. It is commonly proposed that the barren to sub-fertile intrusions are diorite and granodiorite [36,40,41], while fertile intrusions are mostly quartz monzonite, tonalite, and minor granodiorite [42][43][44]. Both granodiorite and granite are considered as the hosts of copper mineralization in the Bondar Hanza ( Figure 3A), in agreement with sub-fertile mineralization.
Four alteration stages are recognized in the Bondar Hanza copper deposit based on the paragenesis and texture relationships, including (i) potassic, (ii) phyllic, (iii) propylitic, and (iv) argillic alteration zones (Figure 3). Potassic alteration is considered by the existence of biotite, quartz, magnetite+albite, and molybdenite, whereas the phyllic zones with quartz, sericite, chlorite, and minor calcite are associated with pyrite and chalcopyrite veinlets ( Figure 3B). The propylitic alteration zones are described by the existence of chlorite, epidote, calcite, pyrite, Fe-oxides, rutile, and albite. The argillic alteration is composed of clay minerals that replace primary aluminosilicate phases, such as K-feldspar, plagioclase, and clay minerals.
In order to identify whether alteration, rock groups, or zones regulate the copper in this region, the correlation coefficient between copper and each of these three categories was calculated. The results indicated that the correlations between the copper and the alteration group, the rock group, and the zone group, as discussed in Table 1, are 7.66%, 39.60%, and 13.05%, respectively. A high correlation coefficient between the copper and the rock group implies that the latter regulates the spatial variability of the former. This is also in agreement with the geological description in this deposit, as discussed above.
Therefore, we decided to take into account the rock groups as the estimation domain for estimating the copper in this deposit.
Further analysis is also necessary to check whether the copper shows different behavior in different rock groups. A boxplot and lognormal probability plot are shown in Figure 7 to check the distribution of copper in rock group A (RG-A), rock group B (RG-B), and rock group C (RG-C). As can be observed, RG-C contains a high amount of copper, contrary to RG-A, which possesses the least amount of copper. This can be seen through the boxplot, where the copper shows the lower median in RG-A compared to the copper in RG-C. The rock groups based on the copper content can be ordered, from most to least important, as RG-C, RG-B, and RG-A. The lognormal probability plot also illustrates that the global distribution of copper in each rock group is different.
Comparison of the global distribution, however, may not be enough to ensure that the behavior of copper is different in each rock group. For more clarification on this, the local distribution of copper was compared in different rock groups as well. Two applicable tools, the madogram and rodogram, are specific measures of spatial variability. The power of h-increments in a traditional variogram is lower in a madogram and rodogram, i.e., 1 and 1 /2 , respectively [25]. These two functions are less sensitive to the extreme values, where they can be used to infer features such as ranges and anisotropy [15]. Therefore, we used these two tools to verify the local continuity of the copper in different rock groups. Figure 8 shows that the copper demonstrates different ranges and local variances in different rock groups. Among them, the copper in RG-A shows a very short-scale range, a medium-scale range is shown in RG-C, and a long-range scale is shown in RG-B. Comparison of the global distribution, however, may not be enough to ensure that the behavior of copper is different in each rock group. For more clarification on this, the local distribution of copper was compared in different rock groups as well. Two applicable tools, the madogram and rodogram, are specific measures of spatial variability. The power of h-increments in a traditional variogram is lower in a madogram and rodogram, i.e., 1 and ½, respectively [25]. These two functions are less sensitive to the extreme values, where they can be used to infer features such as ranges and anisotropy [15]. Therefore, we used these two tools to verify the local continuity of the copper in different rock groups. Figure 8 shows that the copper demonstrates different ranges and local variances in different rock groups. Among them, the copper in RG-A shows a very short-scale range, a medium-scale range is shown in RG-C, and a long-range scale is shown in RG-B.   Comparison of the global distribution, however, may not be enough to ensure that the behavior of copper is different in each rock group. For more clarification on this, the local distribution of copper was compared in different rock groups as well. Two applicable tools, the madogram and rodogram, are specific measures of spatial variability. The power of h-increments in a traditional variogram is lower in a madogram and rodogram, i.e., 1 and ½, respectively [25]. These two functions are less sensitive to the extreme values, where they can be used to infer features such as ranges and anisotropy [15]. Therefore, we used these two tools to verify the local continuity of the copper in different rock groups. Figure 8 shows that the copper demonstrates different ranges and local variances in different rock groups. Among them, the copper in RG-A shows a very short-scale range, a medium-scale range is shown in RG-C, and a long-range scale is shown in RG-B.

Contact Analysis
The previous analysis, to investigate either the local or spatial variability of copper, shows that the underlying continuous variable in this study express different characteristics in different rock groups. Therefore, as a common practice in traditional mineral resource modeling, following a hierarchical modeling technique, the rock groups should be first modeled, and the copper needs to be modeled independently in each built rock group using the data that belongs to the rock group sought. However, this method entails that the behavior of the continuous variable across the boundary of rock groups is hard. In order to investigate this characteristic, a contact analysis [45] is implemented. This tool deals with arrangement of data belonging to one rock type and considers its distance from the border with another rock type. To do so, one usually considers the pairs of data points such that the tail is in one rock type and head is in another rock type. The contact analysis can then be calculated by plotting the mean grade of each rock group as a function of distance to the border between two rock types [3]. Accordingly, the average copper variation in different rock groups was examined versus the distance to their boundaries. The transition of copper was soft (smooth) when crossing the RG-B and RG-C through RG-A, and was hard (abrupt) when crossing RG-C through RG-B; therefore, there is a mixture of hard and soft boundaries in this deposit (Figure 9). type. The contact analysis can then be calculated by plotting the mean grade of each rock group as a function of distance to the border between two rock types [3]. Accordingly, the average copper variation in different rock groups was examined versus the distance to their boundaries. The transition of copper was soft (smooth) when crossing the RG-B and RG-C through RG-A, and was hard (abrupt) when crossing RG-C through RG-B; therefore, there is a mixture of hard and soft boundaries in this deposit ( Figure 9).

Results
Since there is a mixture of hard and soft boundaries in the spatial variation of copper in this deposit, using hierarchical/cascade block modeling is not encouraged, since this method assumes that the transition of continuous variables when crossing the boundary of categorical variable is perfectly hard, leading to production of a 3D spatial model of the deposit with this characteristic. From another point of view, we earlier discussed that the local and spatial behavior of copper in each rock group is different. To solve this problem, following the algorithm presented previously for these particular cases, the mixed (indicator/continuous) cokriging paradigm was used to model the copper and rock groups simultaneously. For this purpose, the rock group, which is a categorical var-

Results
Since there is a mixture of hard and soft boundaries in the spatial variation of copper in this deposit, using hierarchical/cascade block modeling is not encouraged, since this method assumes that the transition of continuous variables when crossing the boundary of categorical variable is perfectly hard, leading to production of a 3D spatial model of the deposit with this characteristic. From another point of view, we earlier discussed that the local and spatial behavior of copper in each rock group is different. To solve this problem, following the algorithm presented previously for these particular cases, the mixed (indicator/continuous) cokriging paradigm was used to model the copper and rock groups simultaneously. For this purpose, the rock group, which is a categorical variable, was converted to indicators. To do so, code 1 was assigned to the location where there is a corresponding rock group; otherwise, 0 was assigned. The same procedure was repeated for all rock groups. As a result, four variables-three indicator variables, i.e., RG-A, RG-B, and RG-C, and one continuous variable (i.e., copper)-were taken into consideration for the practice of modeling. Establishing a cokriging system needs the cross-variogram or cross-covariance to be inferred. For this purpose, first, the variogram of all variables was examined in different directions to identify the potential anisotropies ( Figure 10). The results indicated that the ranges and sills are almost identical in different directions; thus, no specific directional continuity can be recognized in the region. iable, was converted to indicators. To do so, code 1 was assigned to the location where there is a corresponding rock group; otherwise, 0 was assigned. The same procedure was repeated for all rock groups. As a result, four variables-three indicator variables, i.e., RG-A, RG-B, and RG-C, and one continuous variable (i.e., copper)-were taken into consideration for the practice of modeling. Establishing a cokriging system needs the cross-variogram or cross-covariance to be inferred. For this purpose, first, the variogram of all variables was examined in different directions to identify the potential anisotropies ( Figure 10). The results indicated that the ranges and sills are almost identical in different directions; thus, no specific directional continuity can be recognized in the region. Therefore, we considered an omnidirectional covariance for the subsequent analysis of each variable since the simple cokriging system of equations is based on covariance measures. Figure 11 shows the experimental direct and cross-covariance, where a two-structured theoretical covariance with nearly a high nugget effect is fitted accordingly based on a linear model of coregionalization [15,26]. This method of inference en- Therefore, we considered an omnidirectional covariance for the subsequent analysis of each variable since the simple cokriging system of equations is based on covariance measures. Figure 11 shows the experimental direct and cross-covariance, where a two-structured theoretical covariance with nearly a high nugget effect is fitted accordingly based on a linear model of coregionalization [15,26]. This method of inference entails that all covariances have the same ranges, albeit different permissible matrices of sills: Once the covariance models are inferred, one can establish the cokriging system. In this study, a simple cokriging system was used for estimating the copper and for three indicators at grid nodes with a dimension of 5 m × 5 m × 5 m and a moving neighborhood with a range of 2000 m for a search ellipsoid of up to 500 conditioning data. In order to compare the results of this study with those of another approach, a second case was considered where the copper and each indicator were modeled separately by only their direct covariances. In this case, no cross-covariance was considered, and no dependency among the variables was therefore taken into account. Accordingly, two cases were compared: Case I: simultaneous modeling by simple cokriging based on inferred direct and cross-covariances of copper and three indicators-the proposed approach.
Case II: independent modeling by simple kriging based on only a direct covariance of each variable-the traditional approach.
In each of these cases, we were faced with an order relation problem that is a common issue in all indicator cokriging paradigms. In order to circumvent this problem, bounds corrections and normalization [26] were implemented as a postprocessing step so that the summation of the conditional probabilities would reach 1. In the end, the most probable rock group was identified at each target node as a typical act in all indicator cokriging mappings. The results are displayed in Figure 12. As can be seen, the categorical map obtained in Case II, where only simple indicator kriging is applied, is patchy and in some areas unstructured. In addition, the copper map obtained from Case II is very smooth, which may be the case when crossing RG-B and RG-C through RG-A, but it is not the case when crossing RG-C through RG-B. In contrast, the results obtained from Case I, where a cokriging system is applied, are more reasonable in terms of better producing the structure of rock groups in the region. In addition, the copper map obtained from this case seems to report more copper and slightly less variance in some parts of the area. In terms of contact analysis, however, it might be hard to distinguish hard or soft boundaries from the maps by visual inspection. was therefore taken into account. Accordingly, two cases were compared: Case I: simultaneous modeling by simple cokriging based on inferred direct and cross-covariances of copper and three indicators-the proposed approach.
Case II: independent modeling by simple kriging based on only a direct covariance of each variable-the traditional approach.  In each of these cases, we were faced with an order relation problem that is a common issue in all indicator cokriging paradigms. In order to circumvent this problem, bounds corrections and normalization [26] were implemented as a postprocessing step so that the summation of the conditional probabilities would reach 1. In the end, the most probable rock group was identified at each target node as a typical act in all indicator cokriging mappings. The results are displayed in Figure 12. As can be seen, the categorical map obtained in Case II, where only simple indicator kriging is applied, is patchy and in some areas unstructured. In addition, the copper map obtained from Case II is very smooth, which may be the case when crossing RG-B and RG-C through RG-A, but it is not the case when crossing RG-C through RG-B. In contrast, the results obtained from Case I, where a cokriging system is applied, are more reasonable in terms of better producing the structure of rock groups in the region. In addition, the copper map obtained from this case seems to report more copper and slightly less variance in some parts of the area. In terms of contact analysis, however, it might be hard to distinguish hard or soft boundaries from the maps by visual inspection.  Figure 13 illustrates the contact analysis calculated over the resulting copper maps for quantitative comparison. The proposed methodology using cokriging (Case I) is perfectly able to produce the original boundary  Figure 13 illustrates the contact analysis calculated over the resulting copper maps for quantitative comparison. The proposed methodology using cokriging (Case I) is perfectly able to produce the original boundary conditions. The soft boundary for copper crossing RG-B and RG-C through RG-A and the hard boundary crossing RG-C though RG-B were perfectly reproduced; in Case II, all boundaries became soft. This is also somewhat recognizable through the resulting maps of copper presented in Figure 12. Furthermore, the average of copper towards the boundary is different in Case II when crossing RG-B through RG-A, signifying the incompatibility of the resulting maps with the original data ( Figure 13).  Figure 13). An important issue of the cokriging results is the smoothing effect. Frequently, extremely high over-and underestimations of original minimum and maximum values of a distribution, respectively, are unfavorable for the employment of cokriging algorithms in resource estimation. Although, a smoothing effect is inevitable in cokriging results, many scholars are seeking methods to alleviate this problem [46,47]. The histogram of the original copper was compared to the estimated copper obtained from Cases I and II. As can be seen, the smoothing effect was much lower in the estimated values of copper in Case I compared to Case II (Figure 14), although in general its presence is inevitable. The mean and variance of the estimated values are also closer to the original in Case I. This means that the incorporation of rock groups as the secondary variable was highly effective in preserving the high values in the final estimation results, which is a significant improvement over Case II, where An important issue of the cokriging results is the smoothing effect. Frequently, extremely high over-and underestimations of original minimum and maximum values of a distribution, respectively, are unfavorable for the employment of cokriging algorithms in resource estimation. Although, a smoothing effect is inevitable in cokriging results, many scholars are seeking methods to alleviate this problem [46,47].
The histogram of the original copper was compared to the estimated copper obtained from Cases I and II. As can be seen, the smoothing effect was much lower in the estimated values of copper in Case I compared to Case II (Figure 14), although in general its presence is inevitable. The mean and variance of the estimated values are also closer to the original in Case I. This means that the incorporation of rock groups as the secondary variable was highly effective in preserving the high values in the final estimation results, which is a significant improvement over Case II, where no geological information was used in the resource meddling of copper. Detailed statistical parameters of estimated copper in each rock group are also reported in Table 3.
Minerals 2021, 11, x FOR PEER REVIEW 18 of 23 no geological information was used in the resource meddling of copper. Detailed statistical parameters of estimated copper in each rock group are also reported in Table 3.

Validation
In order to check the conditional unbiasedness of the proposed algorithm, a cross-validation [15] was implemented as well. In this approach, a sample is removed and by using other samples, the location of the removed sample was estimated. Taking into account the covariance parameters in both Case I and Case II, the values at the location of that removed sample were estimated, and this was implemented over all sample locations. The estimated values were then compared with the original values. For copper, we calculated the errors and the scatterplot between estimated and true values. As can be seen from Figure 15, both approaches were conditionally unbiased, and Case I produced a slightly lower mean relative error.

Validation
In order to check the conditional unbiasedness of the proposed algorithm, a crossvalidation [15] was implemented as well. In this approach, a sample is removed and by using other samples, the location of the removed sample was estimated. Taking into account the covariance parameters in both Case I and Case II, the values at the location of that removed sample were estimated, and this was implemented over all sample locations. The estimated values were then compared with the original values. For copper, we calculated the errors and the scatterplot between estimated and true values. As can be seen from Figure 15, both approaches were conditionally unbiased, and Case I produced a slightly lower mean relative error.
It was important to determine whether the proposed algorithm is unbiased when predicting copper. Another examination was performed over the rock groups through a confusion matrix. This matrix, also known as an error matrix, is a specific table that enables one to validate the functioning of any predictor algorithm. Basically, it is a contingency table, with two dimensions (i.e., "actual" versus "predicted") that allows one to identify the correct classifications and misclassifications of the categories. In this study, we also used this matrix to explore how many times the proposed algorithm misclassified the actual rock groups at each sample location to assess the accuracy of the proposed approach for assessing its ability to correctly estimate the rock groups. Since true rock groups at sample locations are also available through cross-validation, one needs to derive the predicted rock types at these locations to establish the confusion matrix. The results are presented in Table 4. The accuracy identifies how often the proposed algorithm is correct in predicting rock groups, and precision predicts how often the proposed algorithm is correct in predicting each rock group individually. They are high in both cases, but slightly better in Case I. It was important to determine whether the proposed algorithm is u biased when predicting copper. Another examination was performed o the rock groups through a confusion matrix. This matrix, also known as error matrix, is a specific table that enables one to validate the function of any predictor algorithm. Basically, it is a contingency table, with t dimensions (i.e., "actual" versus "predicted") that allows one to ident the correct classifications and misclassifications of the categories. In t study, we also used this matrix to explore how many times the propos algorithm misclassified the actual rock groups at each sample location assess the accuracy of the proposed approach for assessing its ability correctly estimate the rock groups. Since true rock groups at sample lo tions are also available through cross-validation, one needs to derive predicted rock types at these locations to establish the confusion matr The results are presented in Table 4. The accuracy identifies how often proposed algorithm is correct in predicting rock groups, and precis predicts how often the proposed algorithm is correct in predicting ea rock group individually. They are high in both cases, but slightly better Case I.

Recovery Functions
An important task in mineral resource modeling is calculating the recovery functions. These parameters are the fraction of tonnage, the metal quantity, and the mean grade above different cut-offs. Cut-off grade is the minimum grade required in order for a mineral to be economically mined. Fraction of tonnage is the proportion of values greater than or equal to cut-off, and it represents the fraction of recoverable material above that specific cut-off. Grade-tonnage curves are useful tools for reporting the mineral inventory information. Fraction of tonnage and mean grade above different cut-offs are two important curves. The third graph, quantity of metal versus cut-off grades, can be obtained from the other two.
This information is useful for short-and long-term mine planning and is a basis for economic considerations and evaluations of economic block value, the first step in attaining an optimum mine plan for a project, where the net present value (NPV) is optimum. For this purpose, these functions are calculated with different cut-offs for the whole deposit. The results indicate that Case I, where the proposed cokriging approach is used, produces much better tonnage, metal quantity, and mean grades, when compared to Case II ( Figure 16). These investigations also show that the results obtained from Case I are closer to the global recovery functions that were obtained from original copper values at drillhole data. This signifies that the proposed approach not only is unbiased, but also, in general, produces higher recovery functions, particularly higher mean grade, for further economic evaluations. For instance, the fraction of tonnage, mean grade, and metal quantity above 0.3% Cu for Case I and Case II are: 8.4% and 3.4%, 0.453% and 0.407%, and 0.038% and 0.014%, respectively.
indicate that Case I, where the proposed cokriging approach is used, produces much better tonnage, metal quantity, and mean grades, when compared to Case II ( Figure 16). These investigations also show that the results obtained from Case I are closer to the global recovery functions that were obtained from original copper values at drillhole data. This signifies that the proposed approach not only is unbiased, but also, in general, produces higher recovery functions, particularly higher mean grade, for further economic evaluations. For instance, the fraction of tonnage, mean grade, and metal quantity above 0.3% Cu for Case I and Case II are: 8.4% and 3.4%, 0.453% and 0.407%, and 0.038% and 0.014%, respectively.

Conclusions
A cokriging approach is presented in this study for resource modeling. The algorithm is examined through a copper porphyry deposit in

Conclusions
A cokriging approach is presented in this study for resource modeling. The algorithm is examined through a copper porphyry deposit in Iran, and the results reveal that the method is capable of producing both hard and soft boundaries in the final model. In addition, the proposed technique is less influenced by the smoothing effect when compared to independent kriging, where each variable is modeled separately. This was also confirmed through an investigation of recovery functions. Furthermore, both techniques were statistically assessed through cross-validation, where the outputs showed that both approaches are conditionally unbiased, although the results of the proposed algorithm based on cokriging produced slightly less error. The resulting maps of the geological domain make more sense, and the problem of patchy and unstructured features are resolved by the mixed continuous/indicator cokriging technique. For implementing the mixed cokriging algorithm and producing the figures, we used in-house MATLAB programs.
However, there is some room for improvement in the algorithm. The resulting estimates of indicators still suffer from the order relation problem, and using a continuous covariate could not effectively tackle this challenge. To do so, using compositional kriging [48] is a possible option. Another avenue of research is to use more categorical variables in the process of modeling. Cokriging can simply be extended to include more geological domains. For instance, the mineralized zone, alternation, and geometallurgical units can also be considered as categorical covariates in establishing the cokriging system. However, their cross-dependencies should be investigated prior to any further action. Further research can also involve ordinary cokriging paradigm [49] for testing the proposed algorithm.