Resource Estimation in Multi-Unit Mineral Deposits Using a Multivariate Matérn Correlation Model: An Application in an Iron Ore Deposit of Nkout, Cameroon

Abstract

Modeling the spatial dependence structure of metal grades in the presence of soft boundaries between geological domains is challenging in any mineral resource estimation strategy. The aim of this work was to propose a structural model adapted to this type of geological boundary, based on a multivariate Matérn model that fits the observed direct (within domain) and cross (between domains) correlation structures of metal grades. The methodology was applied to a case study of an iron deposit located in southern Cameroon. Cross-validation scores show that accounting for the grade correlation across domain boundaries improved the traditional workflow, where the grade was estimated in each domain separately. The scores were significantly better when we also ensured that the mean grade was locally invariant from one domain to another to reflect the grade continuity across the domain boundary.

1. Introduction

The partitioning of a mineral deposit into several geological domains makes it possible to incorporate the control that geological information (lithology, mineralization, alteration, and structures, among others) exerts on quantitative variables of interest [1,2,3,4,5,6,7,8,9]. Hereinafter, the use of geological domains will focus on mineral resource (specifically, metal grade) estimation, although domaining is also useful for geotechnical and geometallurgical modeling to estimate variables such as rock quality designation, rock mass rating, work index, or metallurgical recoveries. Geological domaining often incorporates a boundary analysis, which results in either ‘soft’ or ‘hard’ boundaries, depending on whether the metal grades vary continuously across the domain boundaries or not [10,11,12,13,14]. Once the nature of the boundary between the different geological domains has been determined, spatial correlation models must be constructed accordingly and used for grade estimation and/or simulations [10,11,15]. Boundary analysis is therefore an essential stage in exploratory data analysis to diagnose the nature of the geological boundaries. Ignoring the presence of a hard boundary can have a major impact on the grade estimation, in particular, ore/waste dilution, over-smoothing, or local biases in the estimated mineral resources [11,16,17]. Furthermore, depending on whether the interpretation of the boundary analysis leads to a soft or a hard boundary, the spatial correlation model can be very different.

In practice, several tools contribute to diagnosing the nature of a geological boundary [13,14,18]. The hard boundary is defined as an abrupt or sudden variation of a quantitative variable (in our case, a metal grade) when moving from one domain to another. The accepted industry workflow then recommends carrying out a geostatistical analysis in each domain independently of the others. In this way, the error variance as well as other cross-validation scores are often improved compared to analyzing the deposit as a whole [5,10,19,20,21,22]. However, in practice, the transition between two domains is not always abrupt, and a metal grade can evolve continuously across a domain boundary, leading to a soft boundary [1,17,18]. In such a case, accounting for the grade continuity or for its spatial correlation across a boundary may considerably improve the precision of the estimation [2,5].

How to optimally perform the estimation of mineral resources when the boundaries between geological domains are soft remains an open problem. An intuitive and constructive way is to use a linear coregionalization model to estimate the grades using data from adjacent geological domains [10]. Thus, data can be taken on either side of a boundary to perform the estimation at non-sampled points in a given geological domain. Other authors [1,2,5,17] propose dilating the geological domain to include samples from adjacent domains up to a given radius from the boundary of the domain of interest. Another way to better capture the behavior of metal grades across soft boundaries is to jointly estimate [5,23,24] or simulate [6,25,26] the grade and the domain indicators representing the geological information. When applied to real case studies, all these approaches have shown to provide better results than those obtained by classic (independent) domain estimations assuming hard boundaries, or those that do not take domains into account, with a mean squared error that could be reduced by 15% to 30% near the domain boundaries [5,10,17,25].

Although there is a measurable grade correlation across domain boundaries, the statistical and structural properties (e.g., the correlation ranges or the variogram sills and shapes) may differ from one domain to another, indicating a non-stationary behavior. From a modeling point of view, one can distinguish two types of non-stationarity within any data set: non-stationarity in terms of the mean (first-order non-stationarity) and non-stationarity in terms of spatial correlation structure (second-order non-stationarity), for which the mean value and the covariance or the variogram vary when moving from one geological domain to another, although they can remain constant (stationary) inside each single geological domain. For example, the non-stationary covariance class proposed in [7] could allow the model to adapt to spatial variations in the correlation structure by setting a different correlation range of the grade in each geological domain, while accounting for correlations between different domains.

An alternative is to view the non-stationary covariance modeling problem as a multivariate problem, where it is of interest to define a direct covariance for the grade within each geological domain and cross-covariances between domains, e.g., through the well-known linear model of coregionalization [10] or the so-called multivariate Matérn model [27,28]. The latter model has proven to be more versatile than the traditional linear coregionalization model, as it allows different correlation ranges, sills, and shapes for the direct (within each domain) and cross (between domains) covariances, while keeping the number of parameters as small as possible [29].

In this context, the aim of this work was to study the impact on mineral resource estimation based on a multivariate Matérn model to account for soft geological boundaries. The outline of this paper is as follows: Section 2 describes the methodology adopted in terms of boundary analysis and modeling. In Section 3, a case study is carried out, and conclusions are presented in Section 4.

2. Methodology

2.1. Boundary Analysis

Before any study intended to estimate, simulate, or characterize in-situ resources of a mining project on the basis of sampling information collected during exploration, it is important to contextualize each quantitative variable of interest in relation to the geological data (lithology, mineralization, alteration, structures, among others). This helps to build the best possible spatial correlation models for future estimations. Putting each variable in its geological context, one obtains a set of sub-databases that can then undergo independent and/or joint analyses [2,30,31,32].

Quite often in mineral resource estimation, the boundary between geological domains is soft. This implies a gradual or progressive transition of the metal grades from one domain to another. It is therefore essential to build coregionalization models capable of handling such a behavior [10,11,33,34,35].

Two characteristics are essential to determine the spatial behavior of a quantitative variable (here, a metal grade) when crossing a geological boundary. They are the variations of its average value (first-order moment) as a function of the distance from the boundary, and the existence of cross-correlations (second-order moment) across the boundary. In this respect, several exploratory and structural analysis tools have been developed in the literature to study these characteristics:

• correlograms [36] makes it possible to measure the spatial cross-correlations between the grades measured in two different domains;
• cross-to-direct variogram ratios [13] measure the variations in the average grade near the boundary of a geological domain;
• pseudo cross-variograms [17] can also measure the difference in the average grade when crossing the boundary between two different domains;
• lagged scatter plots [14] compare grades taken in two different domains separated by a given distance, allowing the identification of both variations in the average grade and cross-correlations across the boundary.

2.2. Structural Analysis and Modeling

Direct and cross-covariance and/or variogram models make it possible, on the one hand, to characterize the structural behavior of the grade in one geological domain and, on the other hand, to evaluate the correlation that exists between grades taken into different domains. Rather than a global stationary model intended to characterize the spatial dependence structure over the entire mineral deposit under study, a non-stationary model would better capture local structural variations and cross-correlations between geological domains.

Consider a mineralization whose behavior varies from one direction to another (anisotropy), and from one point to another as well (non-stationarity). It would be difficult to infer such a behavior from a set of sparse sampling data, and therefore to construct a global covariance model reflecting this behavior. Therefore, the idea of partitioning the field of study into n domains on the basis of geology offers the possibility of issuing the hypothesis of at least local stationarity within each domain and non-stationarity across domains. Based on this, it is possible to formulate a globally non-stationary model. A way to build such a global model is to adapt the model given in [37] that proposes a non-stationary covariance expression for which the scale factor can vary in space (in the present case, one would restrict this general model so that the scale factor is constant in each geological domain, but different between one domain and another). Applied to a Matérn covariance function with variance σ² and shape factor ν [38], this provides a model where the covariance between the grade measured at location x belonging to domain i (associated with a scale factor α_i > 0) and the grade measured at location x + h belonging to domain j (associated with a scale factor α_j > 0) is of the form

$$\forall \mathit{h}\in {ℝ}^3,\forall i,j\in \left\{1,\dots ,n\right\}, C_{ij}\left(\mathit{h}\right)=M\left(\mathit{h};{\alpha }_{ij},\nu ,{\sigma }_{ij}\right)$$

(1)

with

$${\alpha }_{ij}{}^{-2}=\frac{{\alpha }_i{}^{-2}+{\alpha }_j{}^{-2}}{2}$$

(2)

$${\sigma }_{ij}={\sigma }^2{\left(\frac{{\alpha }_{ij}{}^2}{{\alpha }_i {\alpha }_j}\right)}^{\frac{3}{2}}$$

(3)

and

$$M\left(\mathit{h};\alpha ,\nu ,\sigma \right)={\sigma }^2\frac{2^{1-\nu }}{\mathsf{\Gamma }\left(\nu \right)}{\left(\alpha \Vert \mathit{h}\Vert \right)}^{\nu }{K}_{\nu }\left(\alpha \Vert \mathit{h}\Vert \right)$$

(4)

where Γ stands for the gamma function, K_ν for the modified Bessel function of the second kind of order ν, and $\Vert ·\Vert$ for the Euclidean norm. However, such a model is restrictive, insofar as the sills of C_ii and C_jj are the same (σ²), while the sill σ_ij and the scale factor α_ij of C_ij are fully characterized by the sill and scale factors of C_ii and C_jj.

Another way to account for direct and cross correlations simultaneously is to use a multivariate Matérn model, as proposed in [27]. The latter is more flexible in the modeling of the cross-structure, insofar as the direct (in each domain) and cross (between domains) covariances can have different scale factors, sills and shapes, without the aforementioned restrictions on the scale factors and sills of the direct and cross-structures. The idea is therefore to define an n × n matrix-valued covariance function C(h), whose (i,j)-th entry C_ij is the direct (if i = j) or cross (if i ≠ j) covariance of the grades between domains i and j, of the following form:

$$\forall \mathit{h}\in {ℝ}^3, \mathit{C}\left(\mathit{h}\right)={\mathit{C}}_0 \delta \left(\mathit{h}\right)+\mathit{M}\left(\mathit{h};\mathsf{\alpha },\mathsf{\nu },{\mathit{C}}_1\right)$$

(5)

where C₀ is a positive semidefinite matrix of size n × n, δ(h) is a pure nugget effect model, and M(h; α, ν, C₁) is a multivariate Matérn covariance with matrix-valued parameters α (scale factors), ν (shape factors), and C₁ (sills):

$$\forall \mathit{h}\in {ℝ}^3, \mathit{M}\left(\mathit{h};\mathsf{\alpha },\mathsf{\nu },{\mathit{C}}_1\right)={\mathit{C}}_1\frac{2^{1-\nu }}{\mathsf{\Gamma }\left(\mathsf{\nu }\right)}{\left(\mathsf{\alpha } \Vert \mathit{h}\Vert \right)}^{\mathsf{\nu }}{K}_{\mathsf{\nu }}\left(\mathsf{\alpha } \Vert \mathit{h}\Vert \right)$$

(6)

with all the matrix operations (power, product, ratio, etc.) being elementwise. A wealth of sufficient conditions on (α,ν, C₁) for this model to be a valid matrix-valued covariance are given in [39].

In particular, taking all the shape factors equal to 0.5 leads to a multivariate exponential model:

$$\forall \mathit{h}\in {ℝ}^3, \mathit{C}\left(\mathit{h}\right)={\mathit{C}}_0 \delta \left(\mathit{h}\right)+{\mathit{C}}_1\exp \left(-\frac{3\Vert \mathit{h}\Vert }{\mathit{a}}\right)$$

(7)

where a = 3/α (practical ranges) and all the matrix products, exponentials, and ratios are taken elementwise. In this way, one obtains a globally non-stationary model where the structural parameters (practical ranges a(i,j), sills C₁(i,j), and nugget effects C₀(i,j)) vary from one geological domain to another.

The isotropic models in Equations (6) and (7) can be further generalized to anisotropic models by a rotation and scaling of the spatial coordinates (geometric anisotropy) [40].

2.3. Mean Value Modeling

As usual in geostatistical modeling, the mean value of the grade is assumed to be unknown, which allows this mean to vary globally while being locally constant (at the scale of the moving neighborhood used for estimation). However, it is important to specify the relationship between the mean value in a given domain and that of another domain, which depends on the nature of the domain boundary:

• no relationship (the mean values are different and unrelated) in case of a hard boundary, to reflect an abrupt change of the average grade when crossing the boundary;
• equality (same mean values) in case of a soft boundary, to reflect the continuity of the grade across the boundary.

The latter option amounts to using a variant of ordinary cokriging, called cokriging with related mean values [41], when estimating the grades in domains with soft boundaries. One novelty is the combination of this cokriging variant with the multivariate Matérn covariance (Equation (6)), which is more flexible and parsimonious than the classical linear model of coregionalization, to describe the spatial correlation structure.

3. Case Study: Nkout Center Iron Ore Deposit

3.1. Presentation of the Deposit and Exploratory Data Analysis

The Nkout area (southern Cameroon) has been the subject of particular interest from mining companies and geoscience researchers because of its economic potential for iron mining. Several scientific works that were carried out there led to the delimitation of a deposit whose iron grades and tonnage justify an economically viable exploitation. Overall, the deposit lays on an Archean to Proterozoic cratonic basement that forms the northern part of the northern extension of the Congo craton and includes the Ntem complex, the Dja series, and the Nyong. The Nkout iron zone is located in the Ntem complex. Recent geological studies showed that this area (Nkout) is an oxide iron formation comprising fresh magnetite banded iron formation (BIF) at depth, which weathers and oxidizes towards the surface-forming caps of enriched hematite/martite–goethite ores [42,43,44,45,46].

The Nkout iron deposit is made up of three prospects, namely Nkout East, Nkout Center, and Nkout West. The interest of the present study relates to Nkout Center, an iron deposit formed of itabirites (oxidized facies of rich magnetite of the BIF type) [45]. The mineralization of Nkout Center consists of cycles of finely laminated, crystalline to massively banded magnetite grading into coarsely crystalline, sheared magnetite aggregates contained within a coarse-grained quartz–feldspath matrix often containing amphibole and garnet. The magnetite ore horizons are bound by strongly foliated quartz–biotite amphibole gneiss and metasediment. The mineral species, compositions, mineral associations, and liberation have been studied using automated mineralogy (QEMSCAN) combined with whole rock geochemistry and mineralogical techniques [42].

The available data for this study consist of 329 regularly spaced vertical diamond drill holes in the Nkout Center area, with an average drill core recovery of 42%. The iron grades have been assayed utilizing an industry-standard QA/QC program and composited after assaying at a length of 3 m. The histogram of iron grade data (Figure 1) shows a bimodal distribution. The first mode (the smallest) represents low grades (0–20%), while the second mode represents high grades. The mean iron grade (48.08%) belongs to the zone represented by the second mode.

The analysis of iron grade data (quantitative assay) in relation with the lithology (visual log) suggests that the dataset can be partitioned according to lithological types. Specifically, eight domains (denoted D1 to D8) corresponding to a unique lithology or to similar lithological types were constructed (

Table 1. Geological domain definition.

Domain	Description	Code
Domain 1	Superficial iron-rich laterite and saprolite	D1
Domain 2	Coarse-banded magnetite BIF	D2
Domain 3	Fine-banded magnetite BIF	D3
Domain 4	Undifferentiated rocks	D4
Domain 5	Oxidized rocks (itabirites, hematite–magnetite BIF)	D5
Domain 6	Granitic intrusions, pegmatite	D6
Domain 7	Basal metasedimentary rocks	D7
Domain 8	Gneiss rocks, amphibolite, schist, quartzite	D8

).

Table 2. Statistics of iron grade per geological domain.

Domain	D1	D2	D3	D4	D5	D6	D7	D8	Total
Number of data	886	3022	2985	91	484	287	1378	1520	10,653
Proportion	8%	28%	28%	1%	5%	3%	13%	14%	100%
Minimum	1.69	0.78	3.8	3.96	3.6	0.38	0.38	0.6	0.38
Maximum	66.11	62.69	71.61	65.10	49.73	58.46	60.43	57.63	71.61
Mean	43.07	33.92	44.16	41.55	31.1	12.13	14.56	7.56	48.08
Median	47.23	36.06	44.13	43.53	34.52	7.60	9.84	3.71	44.5
St. deviation	13.90	8.42	13.34	14.12	9.23	12.78	12.45	8.59	10.51

summarizes the main statistical characteristics of the iron grade in each domain and globally. The last three domains (D6, D7, and D8) are ignored in the rest of the analysis because they are waste domains (with a mean iron grade below 15%), and the focus is given to D1 to D5 corresponding to economic ore. Figure 2 presents a plan view of the spatial distribution of these lithological domains and the associated iron grades.

Lagged scatter plots, as well as plots of the average grade as a function of the distance to the geological domain boundaries, were used to perform boundary analysis of the iron grade data in the ore domains (D1 to D5). Their analysis (Figure 3 and Figure 4) suggested that the boundaries between D1–D2, D1–D3, D1–D4, D2–D3, D2–D4, D2–D5 and D3–D4 are soft, with an average iron grade that does not vary or varies little across the boundaries, while D5 has a hard boundary (strong changes in the average iron grade) with D1 and, to a lesser extent, D3 and D4. Overall, D5 seems to be a separate domain and we can construct a group D1–D2–D3–D4 that will be considered jointly and D5 which will be the subject of an individual analysis. This decision will be corroborated during the structural analysis stage (next section), which will reveal that the spatial correlation of the iron grade in D5 is much poorer than in D1 to D4.

Finally, it appears that the data from the central Nkout iron deposit present a combination of several types of boundaries. This leads to two modeling approaches:

• Separate local models: the aim is to characterize the spatial variability within each lithological domain.
• Use a global model to describe the joint behavior of the group D1–D2–D3–D4 with the approach described in Section 2.2 and Section 2.3. Such a model would make it possible to account for variations in the spatial correlation structure of the iron grade when moving from one domain to another.

3.2. Geostatistical Modeling

To study the spatial dependence structure of iron grades, the nonergodic experimental direct (within each domain) and cross (between domains) covariances [47] were calculated along the horizontal and vertical directions, recognized as the main anisotropy directions, for a horizontal lag multiple of 95 m, a vertical lag multiple of 15 m, and a lag tolerance of half the lag. At first sight, the nuggets, sills, and correlation ranges differ with the domain under consideration, with D2 having smaller ranges than the other domains (Figure 5). A nugget + exponential covariance model with a geometric anisotropy was used for fitting each direct or cross covariance for the iron grades in D1 to D4, leading to the parameters of a multivariate exponential model (Equation (7)) indicated in

Table 3. Parameters of fitted multivariate exponential model (D1 to D4).

Index of First Domain (i)	Index of Second Domain (j)	Nugget Effect C0(i,j)	Sill C1(i,j) of Exponential Structure	Horizontal Range (m) of Exponential Structure	Vertical Range (m) of Exponential Structure
1	1	85	108.5	500	166.7
1	2	0	36.88	500	166.7
1	3	0	94.45	500	166.7
1	4	0	32.13	500	166.7
2	2	0	79.33	300	100
2	3	0	46.98	500	166.7
2	4	0	27.22	500	166.7
3	3	43	135.68	500	166.7
3	4	0	57.93	500	166.7
4	4	25	172.70	500	166.7

.

The matrix C₀ is positive semidefinite, which ensures the mathematical validity of the nugget effect component. As for the exponential component, it is a multivariate Matérn model (Equation (6)) with the shape factors matrix (ν) equal to 0.5 times an all-ones matrix and complies with the sufficient condition given in [39] (Theorem 3A):

• ν is conditionally negative semidefinite;
• ν a² (elementwise product and power) is conditionally negative semidefinite;
• C₁a³ν^ν+3/2 exp(−ν)/Γ(ν) (elementwise products, powers, exponentials, and ratios) is positive semidefinite.

Recall that a real symmetric matrix A of size n × n is positive semidefinite if its eigenvalues are nonnegative, and is conditionally negative semidefinite if the symmetric matrix B with entries B(i,j) = A(i,n) + A(n,j) − A(i,j) − A(n,n) is positive semidefinite [48].

Concerning the last domain (D5), a nugget + exponential covariance model with a geometric anisotropy is also used for fitting the experimental covariance, leading to the model parameters indicated in

Table 4. Parameters of fitted univariate exponential model (D5).

Nugget Effect C0	Sill C1 of Exponential Structure	Horizontal Range (m) of Exponential Structure	Vertical Range (m) of Exponential Structure
23.3	61.9	50	70

.

3.3. Cross-Validation

Leave-one-out cross-validation was performed to assess the adequacy of the proposed model to the available data for D1–D4. Three approaches are analyzed comparatively:

• The first approach (hereinafter, approach 1) is a cokriging using the multivariate Matérn model and the additional conditions that the mean values in D1, D2, D3, and D4 are locally the same (but unknown) due to the soft domain boundaries, as explained in Section 2.3. This amounts to a model with stationarity in the first-order moment (the mean value) and non-stationarity in the second-order moment (the covariance, which varies with the domain).
• The second approach (approach 2) uses the same covariance model but gets rid of the previous restriction on the mean values; it is a classical ordinary cokriging where the means are unknown and have no relationships between each other.
• The last approach (approach 3) consists of ordinary kriging performed in each domain separately, i.e., it also gets rid of the cross-covariances as if the domain boundaries were hard.

For each approach, a moving neighborhood of half-size 500 m along the horizontal directions and 200 m along the vertical direction, containing up to 40 data points, was used. When estimating the iron grade of a given domain, classical ordinary cokriging (approach 2) needs data of all the domains for the left-hand side matrix of the cokriging system to be non-singular; however, this inconvenience can be bypassed by removing the grade variables for which no data are available. For instance, when estimating the grade of D1, if no data of D4 are found in the neighborhood, then it suffices to cokrige the grades of D1–D3 (trivariate cokriging), instead of that of D1–D4 (quadrivariate cokriging). There is no such inconvenience with ordinary kriging (approach 3) or with cokriging with related means (approach 1), which only require the moving neighborhood to contain at least one data point to yield a non-singular cokriging system.

Figure 6 shows examples of scatter plots between true and estimated iron grades for approach 1, while

Table 5. Cross-validation scores for the estimation errors in D1 to D4 (full data set).

Approach	Statistics	Domain 1	Domain 2	Domain 3	Domain 4
Approach 1	Number of data	886	3022	2985	91
	Mean	0.192	−0.147	0.000	−0.952
	MAE	5.733	3.787	3.917	6.976
	RMSE	7.691	5.984	5.762	10.017
Approach 2	Number of data	886	3022	2985	91
	Mean error	−0.057	−0.106	−0.012	3.940
	MAE	5.838	3.808	3.934	10.301
	RMSE	8.155	6.067	5.819	15.980
Approach 3	Number of data	886	3022	2985	91
	Mean	−0.019	−0.106	−0.011	3.828
	MAE	5.848	3.810	3.940	10.332
	RMSE	8.192	6.070	5.828	16.117

and

Table 6. Cross-validation scores for the estimation errors in D1 to D4 (subset of data distant no more than 6 m to a geological boundary).

Approach	Statistics	Domain 1	Domain 2	Domain 3	Domain 4
Approach 1	Number of data	282	264	563	33
	Mean	0.920	−0.007	0.086	−1.243
	MAE	4.496	3.477	4.505	6.754
	RMSE	5.909	5.849	6.444	8.970
Approach 2	Number of data	282	264	563	33
	Mean error	0.724	0.386	0.109	9.195
	MAE	4.804	3.719	4.589	15.211
	RMSE	7.526	6.772	6.675	21.523
Approach 3	Number of data	282	264	563	33
	Mean	0.832	0.390	0.124	8.987
	MAE	4.831	3.744	4.615	15.336
	RMSE	7.617	6.815	6.712	21.893

present a summary of the cross-validation scores for all three approaches. Overall, the scores differ from one approach to another, the differences being stronger for the subset of data located close to a domain boundary (Table 6) than for the full set of data (Table 5), as was expected. Approach 1 turns out to be the most accurate. When ignoring the equality relationship between the mean values in the neighborhood of the domain boundaries, the scores significantly deteriorate in D1, D2, and D4 and slightly deteriorate in D3, as observed when comparing approaches 1 and 2: in the subset of data close to a domain boundary (Table 6), the mean absolute error (MAE) increases 7% in D1 and D2, 2% in D3, and 125% in D4, while the root mean squared error (RMSE) increases 27%, 16%, 4%, and 140%, respectively. When also ignoring the cross-correlations of the iron grades across the domain boundaries, there is an extra loss of accuracy, although of small magnitude (between 1% and 2% increase in RMSE for D1 and D4), as seen when comparing approaches 2 and 3.

3.4. Block Models

Block model estimates of iron grade were constructed for each of the approaches under consideration based on the available drill hole data and on a layout of the lithological domains interpreted by mining geologists (Figure 7a). The same cokriging neighborhood as in Section 3.3 was used. Each of the block models is a three-dimensional representation of the spatial distribution of the iron grade in the study area (Figure 7b–d).

Although the theoretical principle underlying the three approaches represented by these figures are different, the block models and the total amount of iron resources are similar (the mean iron grades are 35.00%, 35.09%, and 35.10%, respectively), the differences being noticeable locally in the vicinity of the domain boundaries: approach 1 leads to a block model with smoother grade transitions than approaches 2 and 3, for which the transitions are sharper and the grade model looks more ‘fragmented’, with ‘hard’ geological boundaries; this difference stems from a better use of the sampling data of adjacent domains when estimating the grade within a given domain.

3.5. Discussion

The presented case study highlights two topics that, in our opinion, deserve some discussion. The first one relates to the use of stationary or non-stationary models to represent a regionalized variable such as a metal grade in a mineral deposit, and the second one relates to the benefits or lack of benefits of using cokriging instead of kriging.

As pointed out by numerous authors [14,40,49,50,51], stationarity is a property of the stochastic model (i.e., of the abstract representation of the regionalized variable by a random field), rather than of the regionalized variable itself, and cannot be refuted in general. That is, stationarity is a model decision rather than a model assumption, a decision that, in practice, considerably eases the statistical inference of the model parameters and is motivated by the spatial homogeneity of the data. In mineral resource estimation, however, data homogeneity is questionable when the data behavior strongly depends on the ore type or the rock type, and geological domaining is a commonly accepted practice to deal with heterogeneities.

Our approach builds on usual domaining practice and weakens the stationarity decision regarding both the first- and second-order moments. On the one hand, with the implementation of cokriging in a moving neighborhood, the mean iron grade is locally constant, but it can vary between one neighborhood and another one. On the other hand, the covariance function is constant within each geological domain, but its range and sill can change between one domain and another one. Said another way, one assumes (or rather, one decides) local stationarity for the mean value and ‘piecewise’ stationarity for the covariance function, which makes the spatial structure identification problem approachable with simple multivariate tools and without utilizing fully non-stationary models. For the latter models, the reader is referred to [35,51,52,53,54] and references therein for an overview of different alternatives; note that most of these assume that the variance (covariance at the origin) is constant in space, while our approach can accommodate the case when the variance changes between one geological domain and another one.

Regarding the second topic of discussion, in view of the somehow weak cross-correlations between grades measured in different geological domains (0.49 between D1 and D3, 0.55 between D2 and D3, and less than 0.30 in the other cases, see Figure 3), one can wonder whether cokriging outperforms kriging and is worth the effort. The cross-validation results presented in Section 3.3 proves that there is a significant improvement in the borders of D1, D2, and D4, and a slight improvement for D3. A close look at Table 6 shows that most of the improvement originates from the related mean values (approaches 1 vs. 2), rather than from the cross-covariances (approaches 2 vs. 3). So, even in the case of weak spatial correlations between geological domains, we advocate for the use of cokriging, as it allows accounting for the continuity of the mean grade across soft geological boundaries; a feature that, regrettably, is omitted too often in mineral resource estimation [55].

4. Concluding Remarks

Accounting for geology, the dataset at Nkout Center deposit in southern Cameroon has been divided into several subsets, each corresponding to a particular lithological type or a combination of similar lithological types. A boundary analysis based on lagged scatter plots suggested the coexistence of hard and soft boundaries between lithological domains, with a preponderance of soft boundaries between four out of the five domains of interest. The direct covariances, designed to measure the spatial dependence structure within each domain, revealed dissimilarities in ranges, sills, and nugget effects across the domains. This motivated the modeling of the direct and cross-covariances to represent the joint grade behavior in the different domains, accounting for the correlations across domain boundaries. It should be noted that the only variable under study is the iron grade, but the partitioning into domains leads to different, purely heterotopic, covariates, each associated with a single domain.

In the presence of geological domaining, the simplest idea, widely adopted by practitioners, is to consider each domain individually for geostatistical estimation. However, when the domain boundaries are soft, the spatial correlation model to be built must integrate the local variations within domains and the cross-correlations between domains. The proposed approach takes advantage of the flexibility of the multivariate Matérn covariance model, allowing the ranges, sills, and nugget effects to vary from one domain to another, while minimizing the number of parameters in the coregionalization model. This yields a locally stationary, but globally non-stationary behavior. Cokriging also allows one to account for the fact that there are almost no variations in mean grade when crossing the boundaries of the domains selected to be part of the defined ore group (D1–D2–D3–D4).

The cross-validation scores show that accounting for the grade correlation across boundaries improves the traditional workflow where the estimation is held in each domain separately. The scores are even better when one also ensures that the mean grade is locally invariant from one domain to another.

In conclusion, accounting for the nature of the boundaries between geological domains is essential when estimating mineral resources, and including an iterative boundary analysis before moving on to the spatial analysis step can have a real impact in all mineral estimation domain workflows. In the case of soft boundaries, the assumptions related to both the moment of order 1 (mean value) and the moments of order 2 (direct and cross-covariances) are decisive.

Although cokriging is optimal in terms of being unbiased and minimizing the variance of the estimation error, it generally oversmoothes the reality. As a result, the block model does not reproduce the behavior of the true grades as described by the defined structural model. It would therefore be useful to continue this work by simulating the grades; one can also advocate for simulating the geological domains to account for the uncertainty in their boundaries.