Geostatistical Methodology to Characterize Volcanogenic Massive and Stockwork Ore Deposits

: The Zambujal ore deposit, Neves-Corvo mine, is a zoned volcanogenic sulﬁde deposit of copper and zinc, with massive ores at the top and stockwork ores at the bottom. Metal grades are strongly zoned by ore types. The main methodology hereby proposed combines geostatistical techniques so that an improved stochastic geological model of the Zambujal deposit encompassing morphology and grades is presented. The model of the morphology was made in two main steps. First, a 3D solid of the boundaries was created and then a 3D grid model of the local sulﬁde proportion was simulated. This latter variable was modeled by using rock speciﬁc gravity as a proxy. After that, a conditional grid model of relative copper grades, i.e., recalculated metal grades assuming only the sulﬁde content, is also simulated in accordance with the morphology. At the end, the new tool, metal tonnage cut-off surface, is proposed, which combines copper grades within massive ores and stockwork ores. To validate the results found, the global tonnages of copper obtained by Ordinary Kriging and the proposed methodology are compared.


Introduction
Geological models and metal grade models are fundamental tools used in the mining industry [1]. Both models are generated and constantly updated by taking into account the geological and geochemical information, gathered during the investigations performed at different stages of a mining project. The information provided by the geological modeling helps to better understand the geological behavior of the study area, the spatial distribution of the ores and metal grades, and to further optimize the mining works and mine planning.
The use of geostatistics for modeling mineral deposits boomed with the vulgarization of drilling techniques, sampling and analytical determinations. Despite the increased computational sophistication, traditional geological mapping plays an important role in modern modeling, because it is a heuristic method of passing information to geological models [2,3].
Traditionally, the construction of a mineral deposit model involves the characterization of the morphology and the contents of the elements with economic interest and penalties. In general, the geostatistical evaluation of ore resources encompasses the main steps [1]: (i) definition of the objectives of the study and database inventory, (ii) subdivision of the volume of interest of the mineral deposit into a few regions according to contact with host rock, different ore types, etc. (morphological model), (iii) zoning analysis of the mining variables according to the regions, (iv) creation of a spatial model of each variable under study, (v) estimation of the values of each variable in the non-sampled locations supported by a regular grid of blocks (model of metal grades), (vi) validation of the estimated model, ensuring the consistency of the model with the available data and the geological MoveTM ® software from Midland Valley (Glasgow, UK) was used for geological modeling and visualization. Geostatistical simulations are performed using a high-performance parallel computing version of the DSS algorithm [22,23].

Geological Background
The Neves-Corvo mine is a world-class copper mine in Europe, with a mining area that covers 13.5 km 2 and an exploration concession area of about 112 km 2 . The mine is located in the Alentejo region, about 220 km southeast of Lisbon. It exploits the richest ore deposit ever found in the IPB and contains more than 300 Mt of massive sulfides, which is now considered a giant volcanogenic massive sulfide [19,[24][25][26][27]. Neves-Corvo encompasses six lenticular polymetallic orebodies, Neves, Corvo, Graça, Zambujal, Lombador, and Semblana, which lie on both flanks of a main anticline [28].
The abnormal geochemistry of the Neves-Corvo is unique compared to the other deposits in the IPB: it possesses uncommonly high copper to zinc ratios (100 Cu/(Cu + Zn) > 50%), which more than doubles the other known deposits. The massive sulfide orebodies usually exhibit a classic lateral and vertical metal zonation and rest on the felsic volcanic rocks. As a general rule, the copper-rich ores occur at the base of the massive sulfide lenses, sometimes overlain by zinc-rich ores. Understanding how exhalative volcanic massive sulfide deposits are formed, the environment in which they were formed, and how the metals are saturated during genesis are critical to the exploration and exploitation of these resources.
The Zambujal orebody, which is the focus of the present study, contains significant grades of copper and zinc. It is estimated to be aerially about 600 m by 550 m with a maximum thickness of 55 m [28]. As said, it encompasses a mixture of two ore types, predominant massive ores (RM) in the top part and stockwork ores (RS) (both stringer and impregnation types) in the bottom part.
The massive sulfides mineralization always appears to be intercalated with the black shales, within locations with important silicified content. The massive lenses are connected and constitute single well-defined stratigraphic horizons in Neves-Corvo. The most predominant sulfur mineral is pyrite but chalcopyrite, sphalerite, and galena also occurs above 5%. Concerning the stockwork mineralization, stringer, or fissure ore, it is composed of anastomosing veinlets networks (feeder channels) and disseminations of sulfides that grade into massive sulfide lenses by coalescence of veins and replacements. Stockworks occur principally in highly silicified and chloritized footwall volcanic rocks. Sulfur minerals within stockwork ores are the same as within the massive ores, and host rock occurs in variable proportions. The transition between massive ores and stockwork ores is quick but without being abrupt, and this issue was decisive for the design of the methodology adopted in this study. Figure 1 above display the SG and Cu grades along the depth of one representative borehole intersecting the entire Zambujal orebody, the gradualness of the variations is clearly visible. Distinction between ore types is a prior classification made by the mine company. volcanic rocks. Sulfur minerals within stockwork ores are the same as within the massive ores, and host rock occurs in variable proportions. The transition between massive ores and stockwork ores is quick but without being abrupt, and this issue was decisive for the design of the methodology adopted in this study. Figure 1 above display the SG and Cu grades along the depth of one representative borehole intersecting the entire Zambujal orebody, the gradualness of the variations is clearly visible. Distinction between ore types is a prior classification made by the mine company.

Methodology
A data-driven methodology was designed taking into account the irregular shape of the Zambujal orebody, the blend between ores and the mining company recommendation for modeling a morphological random variable proportion of sulfides within the rock matrix P(x). Stockwork ores exhibit higher recovery of metals in mineral processing as compared to massive ores with equivalent grades, because flotation separation between chalcopyrite and pyrite is lower in efficiency than between chalcopyrite and host rock [21]. Despite the effective copper grades being much lower in the stockwork ores than in the massive ores (internal dilution due to blend of small veins with host rock), the copper grades in the fraction of sulfides within stockwork ores are in average more than the double of those in the massive sulfides, and this issue increases the relevance of the stockwork ores even in regions with poor copper grades.
Therefore, a high-resolution morphological model was built to quantify the local proportion of sulfides for each block of the grid model (the random variable P(x)), followed by a model of relative copper grades (the random variable Y Cu (x)) [29,30]. Both P(x) and Y Cu (x) are regionalized variables and can be modeled by using geostatistical tools [10,11], and the total amounts of metals and ores can be calculated directly from the estimated or simulated values of these variables.
For a sample location x i with a proportion of sulfide minerals P(x i ) and effective grade of copper Z Cu (x i ), the relative grade of copper Y Cu (x i ) is evaluated by (1).
For instance, if a sample had a grade of 2% copper and a P(x) = 0.7, then the relative grade of copper would be equal to 2%/0.7 = 2.86%. By doing so, the metal contents are recalculated and modeled for the sulfide fractions rather than modeling the total copper content.
Modeling of the variable P(x) is no doubt important however it is not measured at borehole cores. Thus, taking into account the principal paragenesis of the deposit, this work proposes the use of the lab-measured SG available for all core boreholes plus the grades of the chemical elements as proxy information for estimating P(x), first at borehole location and then for the entire deposit.
Combining grid models of both variables (P(x) and Y Cu (x)), a parametric resource of copper surface relating the cut-off grade of copper and the proportion of sulfides (morphological constraint) are readily constructed. The ability to construct such a surface in a straightforward way is an important objective and outcome of this work. The proposed methodology encompasses the following main steps:

1.
Build a low-resolution 3D morphological model, delimiting the external boundaries and the main transition between massive ores in the top (RM) and stockwork ores in the bottom (RS). For the construction of this solid model, a classical procedure is used which involves digitalization of closed polygons in several almost parallel cross-sections, interpolation of surfaces, merger of surfaces in order to build a closed volume, and conversion of the volume into a voxel model (grid mining blocks).

2.
Taking into account the paragenesis of the Zambujal deposit, evaluate empirical cdf (s) of the random variable P(x) at boreholes locations conditional to the SG, grades of metals and sulfur grade all obtained from lab measurements to boreholes samples.
To compute local empirical cdf (s) of P(x), and as there is no linear relationship between P(x) and the set of measured variables, and solutions are not unique, we use an approach based on multivariate analysis. First, an extensive table combining all plausible mixtures of main minerals and host rock of the Zambujal paragenesis are made. According to XRD and microscope polished surface analysis already carried out by the mining company to several samples of ores [25], which are in line with the chemical analyzes, only four minerals are present in quantities above 5% that justify the classification proposed in this step: pyrite (FeS 2 , SG = 5. For the sake of illustration; for a mixture of 30% of host rock with, 30% of pyrite and 40% of chalcopyrite, P(x) = 0.70, theoretical computed values are: Fe = 26.14%, S = 30.01%, Cu = 13.85%, and SG = 4.10. Another example, a mixture of 70% of host rock, with 8% of pyrite, 12% of chalcopyrite, 3% of sphalerite and 7% of galena, P(x) = 0.30, will have the following results: Fe = 7.38%, S = 10.39%, Cu = 4.16%, Zn = 2.01%, Pb = 6.06%, and SG = 3.52.
After this extensive theoretical table of mixtures is completed, each set of observed grades of Cu, Zn, Fe, Hg, and S and SG of each core sample of Zambujal are then compared with each row of the theoretical table, and the set of mixtures with the most similar grades and SG are selected. The set of selected rows are able to build a local empirical cdf of the variable P(x). To evaluate the similitude between theoretical values and lab-measured values, due to different variable dispersions, the Gower distance [31] is used.

3.
For estimation plus simulation of n s1 scenarios of this variable, generate n s1 scenarios of P(x) for each borehole sample and compute their average. The Probability Field Simulation (PFS) algorithm [32] is used to draw values from the local cdf (s), instead of Monte Carlo simulation, as it accomplishes with both the local cdf (s) and the spatial continuity model of P(x).

4.
Compute the auxiliary variable relative copper grade Y Cu (x) at sample locations. By using the values generated in step 3, n s1 scenarios of Y Cu (x) plus the average scenario are obtained. 5.
Estimate regional cdf (s) of the studied random variables P(x), Y Cu (x) and SG(x) for both regions RM and RS by Indicator Kriging (IK) and estimate average images of those variables by Ordinary Kriging (OK) [33]. 6.
Build a high-resolution morphological model of P(x) only within the blocks of each region RM or RS, with Direct Sequential Simulation (DSS) conditional to the cdf (s) per ore type region [34]. 7.
Similarly to P(x), simulate the relative grade of Cu Y Cu (x) by using DSS conditional to cdf (s) per ore type region [34]; n s = n s1 ·n s2 realizations are simulated, n s2 for each one of n s1 realizations of P(x).

8.
Calculate ore tonnages, Cu metal tonnages, and average grades with OK results [1,35]. Consider n b the total number of mining blocks, v b the unitary block volume, Z * Cu (x) OK estimates of copper grade, and D * (x) OK estimates of SG, the average grade and total copper tonnage is given by For relative grades, where Y * Cu (x) is the OK estimate of relative copper grade, the same values are given by For simulated relative grades, the calculations are as follows: For sake of validation, both results must be of the same magnitude: T Cu 1 ∼ T Cu 2 ∼ T Cu 3 .

9.
Evaluate uncertainty, computing local variance of the simulated values for P(x) and Y Cu (x). 10. Build a parametric global surface of copper metal tonnages and function of cut-off copper grades in massive ores and in stockwork ores.
Two notes: (1) Variables P(x), Z Cu (x) and Y Cu (x) are not stationary within the Zambujal deposit but are rather constrained to ore types. For this reason, it was essential to delimit RM and RS bodies with 3D computer solids and to constrain the simulations to the regional cdf (s) of the studied variables per RM or RS ore type. A modified version of DSS in which the simulated values are resampled from the regional cdf (s) is used to simulate the entire deposit [34] (steps 6 and 7 of the methodology). (2) The boreholes assays are made mostly at one meter length and the block size is 2 m × 2 m × 2 m (this is the size of the mining blocks used by the mining company, since it is considered the most suitable for stope design and mine planning). In order to make these scales compatible, prior to the application of the OK and DSS, the values of the boreholes are transferred to the blocks of the model that contained them, and when there are several samples in each block the arithmetic mean is made.

Case Study and Results
The case study follows the previous mentioned steps, combining geological heuristic mapping, kriging estimation and stochastic simulation algorithms. The database provided by Lundin Mining compiles information from more than 300 boreholes, for a total of about 58 km of drills. This database contains four tables of data concerning borehole information, survey, geology/ores and assays. With fewer exceptions, the boreholes assays are reported to one meter length intervals. It has implemented a quality assurance/quality control (QAQC) program to monitor sampling selection, sample preparation, and analytical process in order to assess the confidence in the stored data.
For sake of confidentiality, quantities are presented in all figures and tables previously multiplied by a constant.

Low Resolution 3D Morphological Model
The two mineralized domains RM and RS were manually digitized along several parallel cross-sections according to the borehole data and expertize heuristics. The distinction between RM and RS in boreholes used the classification already proposed by the mining company. Then, three surfaces were interpolated: the top of the massive ore (RM), the bottom of the massive ore/top of the stockwork ore (RS), and the bottom of the RS (Figure 2a). The digitalization ensures a closed volume between surfaces. Then, this model of surfaces was converted into a mining block model (block size 2 m × 2 m × 2 m, total number of blocks is 1,505,350) coded by RM or RS (Figure 2b).

Low Resolution 3D Morphological Model
The two mineralized domains RM and RS were manually digitized along several parallel cross-sections according to the borehole data and expertize heuristics. The distinction between RM and RS in boreholes used the classification already proposed by the mining company. Then, three surfaces were interpolated: the top of the massive ore (RM), the bottom of the massive ore/top of the stockwork ore (RS), and the bottom of the RS (Figure 2(a)). The digitalization ensures a closed volume between surfaces. Then, this model of surfaces was converted into a mining block model (block size 2 m × 2 m × 2 m, total number of blocks is 1,505,350) coded by RM or RS (Figure 2(b)).

Evaluate the Empirical Cdf of the Random Variable P(x) at Boreholes Locations
Before the geostatistical steps, the variable P(x) was estimated for all samples following the procedure mentioned in step 2 of Section 2.2. To make this estimate at the sample locations, the measured levels of the most abundant chemical elements (Fe, Cu, Zn, Sn, Pb, and S) were used, as well as the SG(x) of the sampled cores and the main minerals of the deposit (pyrite, chalcopyrite, galena, blend) plus host rock. Table 1 lists the main minerals of Zambujal, chemical formulas, specific gravity, and theoretical chemical element percentages according to the atomic numbers of the elements. It is important to emphasize the difference between the specific gravity of the host rock (about 2.8) and the specific gravity of the sulfide minerals (all above 4). An extensive table listing all admissible combinations of these minerals and the host rock was then made with a resolution of 1% for proportion, and, for each combination, the theoretical values of the element grades (based on the chemical formulas) and the composite specific gravity were computed. Table 2 shows eight examples of the 4,598,126 identified as plausible. Then, the laboratory-measured values of grades and SG are compared with each item of the extensive table, and a set of the closest combination of minerals plus host rock was extracted. Based on a sensitivity analysis for a selected set of 5% of the available assays [28], 10 theoretical mixtures of minerals plus host rock were considered sufficient to characterize the local empirical cdf. Each selected set of 10 values constitutes the local empirical cdf of P(x) variable at the locations of the boreholes. Figure 3a shows the histogram of the averages of the values of P(x) and Figure 3b shows the histogram of the coefficients of variation (dispersion measurement of the cdf ). Values of P(x) clearly show two populations, respectively, for RM and RS ores, and the separation threshold is about 0.6 (up to 60% of sulfides is stockwork, above 60% is massive ores). For the majority of the cdf(s), the coefficient of variance is low (for 80% of the cdf(s) it is below 0.15). The maximum observed range of P(x) (max-min) was 0.1 and the average range of all solutions was 0.04. These values substantiate the proposed estimation of the local P(x) cdf(s).

Estimation and Simulation of P(x), ZCu(x)and YCu(x) Variables
Cdf(s) of P(x) were used for variography, estimation and simulation of this variable within the entire deposit of Zambujal, as well as the calculation of relative grades YCu(x) at borehole locations.
To identify the principal structural directions of the deposit, experimental variograms were computed for more than 20 directions according to the direction of the boreholes for all modeled variables [28]. Besides, for nugget effect evaluation, downhole variograms were also computed.  Values of P(x) clearly show two populations, respectively, for RM and RS ores, and the separation threshold is about 0.6 (up to 60% of sulfides is stockwork, above 60% is massive ores). For the majority of the cdf (s), the coefficient of variance is low (for 80% of the cdf (s) it is below 0.15). The maximum observed range of P(x) (max-min) was 0.1 and the average range of all solutions was 0.04. These values substantiate the proposed estimation of the local P(x) cdf (s).

Estimation and Simulation of P(x), Z Cu (x)and Y Cu (x) Variables
Cdf (s) of P(x) were used for variography, estimation and simulation of this variable within the entire deposit of Zambujal, as well as the calculation of relative grades Y Cu (x) at borehole locations.
To identify the principal structural directions of the deposit, experimental variograms were computed for more than 20 directions according to the direction of the boreholes for all modeled variables [28]. Besides, for nugget effect evaluation, downhole variograms were also computed.  (1) and (3) are contained in the plane of the sections used in the previous step of digitalization, while the direction (3) intersects those sections. The nugget effect displayed by downhole variograms is about 5% for P(x) and 10% for Z Cu (x) and Y Cu (x); because these amounts themselves are relatively low and with the upscaling between the borehole scale and the block model scale, no nugget effect was considered for OK or DSS.
Experimental variograms adjustments are of good quality for all variables. All variograms were fitted with spherical functions (SPH) and no nugget effect was added. Anisotropy ratio is remarkable for variables P(x) and SG(x) about 6×, major range is 300 m and minor range is 60 m. For grades and relative grades anisotropy ratio is much smaller about 1.33×, major range is about 40 m and minor range is 30 m. Figure 4 shows the experimental variograms for P(x) in the above-mentioned directions and Figure 5 the experimental variograms for Y Cu (x) for the same directions.

Estimation and Simulation of P(x), ZCu(x)and YCu(x) Variables
Cdf(s) of P(x) were used for variography, estimation and simulation of this variable within the entire deposit of Zambujal, as well as the calculation of relative grades YCu(x) at borehole locations.
To identify the principal structural directions of the deposit, experimental variograms were computed for more than 20 directions according to the direction of the boreholes for all modeled variables [28]. Besides, for nugget effect evaluation, downhole variograms were also computed. Based on the results, three orthogonal directions: (1) W-E and 20° dip (90°; −20°), (2) N-S, horizontal (0; 0), and (3) W-E with 70° dip (90°; 70°) were selected as principal directions. The directions (1) and (3) are contained in the plane of the sections used in the previous step of digitalization, while the direction (3) intersects those sections. The nugget effect displayed by downhole variograms is about 5% for P(x) and 10% for ZCu(x) and YCu(x); because these amounts themselves are relatively low and with the upscaling between the borehole scale and the block model scale, no nugget effect was considered for OK or DSS.
Experimental variograms adjustments are of good quality for all variables. All variograms were fitted with spherical functions (SPH) and no nugget effect was added. Anisotropy ratio is remarkable for variables P(x) and SG(x) about 6×, major range is 300 m and minor range is 60 m. For grades and relative grades anisotropy ratio is much smaller about 1.33×, major range is about 40 m and minor range is 30 m. Figure 4 shows the experimental variograms for P(x) in the above-mentioned directions and Figure 5 the experimental variograms for YCu(x) for the same directions.   As explained before, the next experimental steps are: (i) draw values from the local cdf (s) of P(x), (ii) estimation of P(x), Z Cu (x) and Y Cu (x) variables by OK, and (iii) DSS of P(x) and Y Cu (x) variables. Thirty realizations of P(x) plus three realizations of Y Cu (x) for each one of P(x) are performed: (n s = 90; n s1 = 30; n s2 = 3). Ninety realizations for Y Cu (x) are common for simulation of variables of this type, and the smooth appearance of the output images of uncertainty confirms this issue. Figure 5 shows a 3D view of a single from the thirty realizations of P(x) at the locations of the boreholes done by PFS. It became perceptible that the higher values are at the top of the structure and the lower values at the bottom.
Thirty realizations of P(x) plus three realizations of YCu(x) for each one of P(x) are performed: (ns = 90; ns1 = 30; ns2 = 3). Ninety realizations for YCu(x) are common for simulation of variables of this type, and the smooth appearance of the output images of uncertainty confirms this issue. Figure 5 shows a 3D view of a single from the thirty realizations of P(x) at the locations of the boreholes done by PFS. It became perceptible that the higher values are at the top of the structure and the lower values at the bottom. The proportion of sulfide ore minerals P(x) was simulated using DSS for the two regions (RM and RS) at once. Instead of using a global cdf, the simulation of P(x) was performed using a modified version of DSS that uses local cdf(s) for each region, P(x)|RM and P(x)|RS. For declustering purposes, P(x)|RM and P(x)|RS are estimated using IK of 20 P(x) classes.
The simulation of relative copper grades YCu(x) follows the same approach as taken for P(x), namely the use of IK to estimate conditional cdf of YCu(x)|RM and YCu(x)|RS, followed by DSS with local histograms.  The proportion of sulfide ore minerals P(x) was simulated using DSS for the two regions (RM and RS) at once. Instead of using a global cdf, the simulation of P(x) was performed using a modified version of DSS that uses local cdf (s) for each region, P(x)|RM and P(x)|RS. For declustering purposes, P(x)|RM and P(x)|RS are estimated using IK of 20 P(x) classes.
The simulation of relative copper grades Y Cu (x) follows the same approach as taken for P(x), namely the use of IK to estimate conditional cdf of Y Cu (x)|RM and Y Cu (x)|RS, followed by DSS with local histograms. From Figures 5, 6(a) and 7b,e,f, higher values of P(x) are located at the top of the structure, and lower values at the bottom, and the transition between the two regions is quick; in fact, histogram of Figure 3 supports this issue. Figures 6(b) and 7c,d,h,i show that the relative grades in stockwork ores are significantly higher than those in the massive ores. The transition of grades between regions is smooth, which is in accordance with the variability of the grades along the boreholes. From Figures 5, 6a and 7b,e,f, higher values of P(x) are located at the top of the structure, and lower values at the bottom, and the transition between the two regions is quick; in fact, histogram of Figure 3 supports this issue. Figures 6b and 7c,d,h,i show that the relative grades in stockwork ores are significantly higher than those in the massive ores. The transition of grades between regions is smooth, which is in accordance with the variability of the grades along the boreholes.
From Figures 5, 6(a) and 7b,e,f, higher values of P(x) are located at the top of the structure, and lower values at the bottom, and the transition between the two regions is quick; in fact, histogram of Figure 3 supports this issue. Figures 6(b) and 7c,d,h,i show that the relative grades in stockwork ores are significantly higher than those in the massive ores. The transition of grades between regions is smooth, which is in accordance with the variability of the grades along the boreholes. For sake of validation, histograms and variograms of selected realizations were computed and the results (Figure 8 show an example of one realization of P(x) and YCu(x)) show that they are quite well reproduced [28]. Also, copper grades were computed from the inverse relationship. * ( ) = * ( ) · ( ) (5) For sake of validation, histograms and variograms of selected realizations were computed and the results (Figure 8 show an example of one realization of P(x) and Y Cu (x)) show that they are quite well reproduced [28]. Also, copper grades were computed from the inverse relationship.
simulated images of P(x); (h) one of the ninety simulated images of YCu(x); (i) average image of YCu(x) calculated from simulations; (j) image of variance (uncertainty) calculated from the set of simulated images of YCu(x); (k) binary image of variance (uncertainty) identifying the 25% of blocks with lower local variance calculated from the set of simulated images of P(x); (l) binary image of variance (uncertainty) identifying the 25% of blocks with lower local variance calculated from the set of simulated images of YCu(x).
For sake of validation, histograms and variograms of selected realizations were computed and the results (Figure 8 show an example of one realization of P(x) and YCu(x)) show that they are quite well reproduced [28]. Also, copper grades were computed from the inverse relationship. * ( ) = * ( ) · ( ) It is important to note that the values obtained with (5) are similar to the estimated ones, and the range is similar. This evidence is a good indicator that the proposed methodology can be applied and does not bias the results.
In addition, for validation of results, the average of grades obtained by using OK and DSS for the massive ore region, the stockwork ore region, and both regions combined is compared in Table 3. Deviations are very small, confirming that the use of relative grades YCu(x) and DSS with local histograms generates unbiased realizations. It is important to note that the values obtained with (5) are similar to the estimated ones, and the range is similar. This evidence is a good indicator that the proposed methodology can be applied and does not bias the results.
In addition, for validation of results, the average of grades obtained by using OK and DSS for the massive ore region, the stockwork ore region, and both regions combined is compared in Table 3. Deviations are very small, confirming that the use of relative grades Y Cu (x) and DSS with local histograms generates unbiased realizations.

Construction of the Tonnage-Surface Function for Copper Grades and Geological Dilution
Using the simulated values of the variables P(x) and Y Cu (x), average estimates for [P(x i )]*, [Y Cu (x i )]*, and [Z Cu (x i )]* = [P(x i )]*·[Y Cu (x i )]* were computed for each mine block x i . The tonnage of copper (T Cu ) within each mine block x i with volume V and ore specific gravity SG(x i ) is given by: A surface representing the total average tonnage of copper as function of copper cut-offs for both massive (RM) and stockwork ores (RS) is now presented in Figure 9a; two vertical planes are also displayed: (blue) equal copper cut-offs for massive and stockwork ores; (green) cut-offs for massive ores are 1.2× of the ones of stockwork ores. The two corresponding cut-off curves, respectively for equal cut-off grades and richer copper cut-offs in RM than is RS (1.2×) are displayed in Figure 9b. The classification of each block as RM or RS ore type was done by applying a threshold value of 60% (0.6) to the morphological model of P(x) (see Figure 3a).
A surface representing the total average tonnage of copper as function of copper cut-offs for both massive (RM) and stockwork ores (RS) is now presented in Figure 9(a); two vertical planes are also displayed: (blue) equal copper cut-offs for massive and stockwork ores; (green) cut-offs for massive ores are 1.2× of the ones of stockwork ores. The two corresponding cut-off curves, respectively for equal cut-off grades and richer copper cut-offs in RM than is RS (1.2×) are displayed in Figure 9(b). The classification of each block as RM or RS ore type was done by applying a threshold value of 60% (0.6) to the morphological model of P(x) (see Figure 3(a)).  The maximum resource is 528,000 tons and is computed for all blocks independently of the copper grade and of being RM or RS ore type. If a Cu cut-off of 1.3% is applied for the entire deposit (RM + RS), the available copper tonnage is 261,905 tons (line in blue at Figure 9b). However, if a Cu cut-off of 1.3% is applied for the RM ore and Cu cut-off of 1.08% is applied for RS ores (line in red at Figure 8b, Z RM Cu (x) = 1.2 Z RS Cu (x)), the available copper tonnage increased to 280,100 tons.

Discussion and Conclusions
The proposed methodology highlights the importance of the stockwork ores in the Zambujal deposit, and as mineral processing is more efficient in stockwork ores, metal cut-offs should take into account the proportion of sulfides relative to host rock. Also, Zambujal stockwork ores present lower grades of penalty elements, such as lead, arsenic, and bismuth.
When building a grade model of a VMS deposit with massive ores and stockwork, it is important to take into account the following issues: 1.
The proportion between sulfides and host rock P(x) can be modeled as a morphological random variable. Local values of P(x) at the borehole locations allow calculating the auxiliary variable relative copper grades Y Cu (x), meaning the hypothetical metal grades if host rock is not considered.

2.
Modeling the relative copper grades variable Y Cu (x) together with the proportion variable P(x), adds information to each mine block and enables the establishment of conditional cut-off grades to grades and morphology.
The variables in this study are modeled by using Ordinary Kriging (OK) and Direct Sequential Simulation (DSS). OK allows estimating the best unbiased total amount of ores and metal for each domain (RM and RS). In the present case study, OK and DSS results are of the same order of magnitude (see Table 3); hence, no bias was made in any calculation.
Given the available information, the realizations generated by the simulation algorithm, allows for an evaluation of the locations with highest uncertainty, taking into account the distance to the closest samples, the spatial continuity and the local heterogeneity of the observations. The Zambujal deposit has sufficient information to show results almost everywhere with confidence.
The proposed methodology is relevant in the sense that geostatistics demonstrated the importance of the stockwork ores in the Neves-Corvo deposits, in particular the Zambujal deposit. The Neves-Corvo deposit is known to possess an uncommonly high copper to zinc ratio. In the deeper parts of the ore masses, corresponding to the stockwork ores, it is important to note the extremely copper rich ores reflect zone refining and late enrichment. One relevant point is the fact that stockwork ores are usually undervalued because high metal grades are diluted by the host rock within this region.