Resource and Reserve Calculation in Seam-Shaped Mineral Deposits ; A New Approach : “ The Pentahedral Method ”

In recent years, the introduction of specific software for the evaluation of geological resources and mineral reserves has increased significantly thanks to the use of block models capable of working with large databases and applicable to virtually all types of deposits. It is only in layered, tabular-shaped deposits where the use of block models poses certain drawbacks, not only in terms of visual representation but also during the process of interpolation. Other calculation methods for tabular bodies such as sections, polygons, and triangles work with 2D projections but not with 3D. The “Pentahedral” method described here is undoubtedly an innovative method that allows work to always be conducted in 3D, providing a fairly accurate representation of tabular bodies and making it possible to carry out calculations of resources and reserves using any interpolation method. It is demonstrated with figures and tables of the Carlés mineral deposit, a well-developed exoskarn gold deposit in northwestern Spain (Asturias) where the authors have worked extensively. The pentahedral method takes into account not only geological and stratigraphic data from the model interpolation, but also mining concepts such as “minimum thickness” related to the minimum seam size that can be economically and technically mined, and “overbreak” related to the dilution effect that appears during the mining process due to over-excavation.


Introduction
Once the geological, geochemical, and structural information of a mineral deposit are known, it must be represented in 3D in order to obtain sufficiently accurate estimates of its mineralogical content [1], which will lead to the estimation of both mineral resources and mineral reserves, as is explained in international standards such as those from the Canadian Institute of Mining, Metallurgy and Petroleum (CIM) [2].
According to CIM standards, a Mineral Reserve is the economically mineable part of a Measured or Indicated Mineral Resource demonstrated by at least a Preliminary Feasibility Study.This study must include adequate information on mining, processing, metallurgical, economic, and other relevant factors that demonstrate, at the time of reporting, that economic extraction can be justified.A Mineral Reserve includes diluting materials and allowances for losses that may occur when the material is mined.
Despite the existence of several resource estimation methodologies [3], from traditional to geostatistical ones [4][5][6] that use the block model as a geometrical representation, in the particular case Minerals 2017, 7, 72 2 of 21 of tabular or seam-shaped deposits these methods do not allow a sufficiently accurate representation or interpolation [7].
Hexahedral block models are widely used to simulate and classify resources.Reference [8] shows a study of 120 recent NI 43-101 mineral deposit technical reports, where kriging was the most elaborate one, although used in only seven of the reports.Other more refined methods, such as the sequential indicator simulation (SIS, a pixel-based model) or indicator kriging, are mainly usable in reservoir simulations.Implicit modeling is starting to be implemented in commercial simulation packages, e.g., Leapfrog Radial Basis Functions, Surpac Dynamic Shells Module or Minesight Implicit Modeler, but is recommended for use as an additional tool to provide initial and quick resource analysis, not as a complete replacement of block models [9].
State-of-the-art software uses algorithms that sequentially divide each block into smaller blocks in order to follow the geological structures (sub-blocking), but this approach does not achieve the level of precision required and results in models with an extremely high number of blocks; see [1].Additionally, sub-blocking uses the composition of the samples surveyed, but does not make an interpolation using the intersections between the surveys and the mineral structure, considering a cutoff, or a minimal thickness, or even considering lateral dilution due to over-excavation.In this paper, it is the "operational mining dilution" that occurs at the time of mining, as defined in [10], which will be referred to as dilution.
Advances have been made in the modeling of local anisotropy in complex geological structures, see [11], but this method does not include mining parameters such as the minimum thickness values or the dilution.
This paper uses a new method to model tabular deposits, as exemplified by the "Zona M" seam of the northern side of the Carlés deposit, in Spain.This deposit is a gold copper skarn, located in the Rio Narcea Gold Belt [12][13][14] formed at the contact between the Carlés granodiorite and the Nieva limestone of the Rañeces Group [15,16].It consists of several seam-shaped mineral bodies, some of them with economical grades of Au and Cu, with thicknesses ranging from 1 to 6 m and dips from 40 • in the north (55 • ) direction (eastern areas) to 70 • in the north (30 • ) direction (northern bodies).
In the Carlés deposit there are several seams, each of them continuous.They are separated by faults or are located in different contact areas.As is standard practice in classical block models, bodies or seams divided by faults should be treated as different bodies and studied separately.
RECMIN (from: "REC ursos MIN eros" which translates as "mineral resources" in English) software will be used throughout the research in order to manage data and visualization; see [17].RecMin is free software used in mineral deposit modeling, developed entirely by the authors over the past 25 years and thoroughly tested in real working conditions at dozens of mining projects.
RecMin code was supplemented with enough algorithms to consider this new calculation method, described further on.This paper will show how to work with calculation units geometrically different from blocks.The pentahedral method also adds information usually not taken into account in block modeling, such as "lateral dilution" or "mining intercept".

Research Objectives
The main goal of this research is the evaluation of the seam-shaped mineralized structure using a new methodology that provides a precise representation of morphology and physical and chemical characteristics, and storing this information in a suitable manner in order to proceed with resource estimations.The specific objectives of the method are: • Three-dimensional (3D) modeling is mandatory; • It must allow different interpolation algorithms, including geostatistics;

•
The calculation detail has to be user-defined; • Dilution (material below the cutoff) should be considered, if necessary; • It must allow estimates with a categorization using distances, number of intercepts, or a combination of both parameters; • It must allow interpolation, taking into account boreholes filtered by distance: neighborhood influence;

•
The results will be stored in a Structured Query Language (SQL) database.

Geological Characteristics and Constraints
The Carlés Cu-Mo-Au ore deposit is located in the Rio Narcea Gold Belt (RNGB), in the Cantabrian zone of the Iberian Massif.It is related to a small postkinematic calcalkaline monzogranite, which intrudes as a cedar tree laccolith into the upper siliciclastic Furada Formation (shales, oolitic ironstones, and, at the top, some carbonate layers of late Silurian age) and the Nieva carbonates (mostly limestones and some dolostones of early Devonian age), developing skarns bearing copper, molybdenum, and gold of economic interest (Figure 1A).The results will be stored in a Structured Query Language (SQL) database.

Geological Characteristics and Constraints
The Carlés Cu-Mo-Au ore deposit is located in the Rio Narcea Gold Belt (RNGB), in the Cantabrian zone of the Iberian Massif.It is related to a small postkinematic calcalkaline monzogranite, which intrudes as a cedar tree laccolith into the upper siliciclastic Furada Formation (shales, oolitic ironstones, and, at the top, some carbonate layers of late Silurian age) and the Nieva carbonates (mostly limestones and some dolostones of early Devonian age), developing skarns bearing copper, molybdenum, and gold of economic interest (Figure 1A).According to the paper by Martin-Izard et al. [16], the Carlés deposit consists mainly of a well-developed exoskarn [18].The Carlés intrusive is a monzogranite, up to 800 m in diameter, which intrudes along the contact of a Silurian siliciclastic formation and a Devonian carbonate formation.It is a light gray, medium-grained, equigranular igneous rock with plagioclase, biotite, quartz, K-feldspar, and magnesiohornblende.The accessory minerals are zircon, rutile, ilmenite, magnetite, and apatite.The plagioclase (An 30-47) has normal zoning with a common narrow albite-oligoclase (An 8-20) rim.The skarn is at the contact between the monzogranite and the carbonate rock, where a Cu-Au-Mo calcic skarn developed.The mineralized area is divided into three main sectors (east, west and north sectors; See Figure 1), all of which have been studied by Martin-Izard et al. [13].At the moment the east and north sectors are being exploited (Figure 1A-C), both by open pit and also by underground mining.
The exoskarn is mostly calcic skarn made up of early garnet and pyroxene, and later amphibole, magnetite, and sulfides (Figure 2A,B).The ore is related to the skarn retrogradation and post-skarn veining and faulting.According to the paper by Martin-Izard et al. [16], the Carlés deposit consists mainly of a welldeveloped exoskarn [18].The Carlés intrusive is a monzogranite, up to 800 m in diameter, which intrudes along the contact of a Silurian siliciclastic formation and a Devonian carbonate formation.It is a light gray, medium-grained, equigranular igneous rock with plagioclase, biotite, quartz, Kfeldspar, and magnesiohornblende.The accessory minerals are zircon, rutile, ilmenite, magnetite, and apatite.The plagioclase (An 30-47) has normal zoning with a common narrow albite-oligoclase (An 8-20) rim.The skarn is at the contact between the monzogranite and the carbonate rock, where a Cu-Au-Mo calcic skarn developed.The mineralized area is divided into three main sectors (east, west and north sectors; See Figure 1), all of which have been studied by Martin-Izard et al. [13].At the moment the east and north sectors are being exploited (Figure 1A-C), both by open pit and also by underground mining.
The exoskarn is mostly calcic skarn made up of early garnet and pyroxene, and later amphibole, magnetite, and sulfides (Figure 2A,B).The ore is related to the skarn retrogradation and post-skarn veining and faulting.The metasomatic phenomena give way to alternating layers, centimeters to meters in thickness, of garnet-dominant and pyroxene-dominant skarn.Retrograde alteration of the garnet-pyroxene skarn to form a pyroxene, amphibole, magnetite and sulfide minerals assemblage is a very important process in the Carlés deposit.Above all, metallic mineralization at Carlés is most important in lenses within amphibole-pyroxene bands.

Methodology
"Zona M" in the northern area of the Carlés deposit is the zone this study will focus on.The calculation procedure was developed in the following phases: 1. Determination of the seam localization from the geological information and survey data.2. Construction of a 3D-modeled surface at the center of the seam.3. Generation of an equally-spaced cloud of points in the seam centered surface.4. Calculation of the intercept values of the surveys crossing the seam. 5. Execution of a geostatistical study of the intercept values.6. Interpolation of the grade and the thickness of each point of the cloud of points, using the intercept information.7. Definition of the estimation categories and generation of the calculation units.8. Database generation and representation.The metasomatic phenomena give way to alternating layers, centimeters to meters in thickness, of garnet-dominant and pyroxene-dominant skarn.Retrograde alteration of the garnet-pyroxene skarn to form a pyroxene, amphibole, magnetite and sulfide minerals assemblage is a very important process in the Carlés deposit.Above all, metallic mineralization at Carlés is most important in lenses within amphibole-pyroxene bands.

Methodology
"Zona M" in the northern area of the Carlés deposit is the zone this study will focus on.The calculation procedure was developed in the following phases: 1.
Determination of the seam localization from the geological information and survey data.

2.
Construction of a 3D-modeled surface at the center of the seam.

3.
Generation of an equally-spaced cloud of points in the seam centered surface.

4.
Calculation of the intercept values of the surveys crossing the seam.

5.
Execution of a geostatistical study of the intercept values.6.
Interpolation of the grade and the thickness of each point of the cloud of points, using the intercept information.
Definition of the estimation categories and generation of the calculation units.8.
Database generation and representation.

Determination of the Seam Localization from the Geological Information and Surveys
For each and every borehole, the intercept interval of the "ZONA M" mineralized body was determined, without considering the grade value, only considering the skarn lithology.
The seam localization has to be broad enough to allow late adjustments as the modeling defines the correct economical limits, from parameters such as cutoff grade, minimal thickness, or lateral dilution.Figure 3 shows graphically the seam limits at a borehole intercept and the associated bar chart displaying grade data.

Determination of the Seam Localization from the Geological Information and Surveys
For each and every borehole, the intercept interval of the "ZONA M" mineralized body was determined, without considering the grade value, only considering the skarn lithology.
The seam localization has to be broad enough to allow late adjustments as the modeling defines the correct economical limits, from parameters such as cutoff grade, minimal thickness, or lateral dilution.Figure 3 shows graphically the seam limits at a borehole intercept and the associated bar chart displaying grade data.

Construction of a 3D Modeled Surface at the Center of the Seam
As stated above, once the borehole intervals that intercept the mineralized body have been defined, a surface at the center of this body has to be created, from where the interpolation point mesh will be created.As noted above, the surface does not need to be located exactly, as the software will adjust its final position according to the real values of the borehole intercepts and the body.
The easiest way to proceed is to work section by section (Figure 4). Figure 4A shows a cross section of boreholes depicted in white whereas the intercepts with the body are depicted in red.The green line connects the centers of each red line section.Figure 4B simultaneously shows several sections and their corresponding centerlines, which are triangulated to generate a surface.Figure 4C shows the final calculated surface, which will henceforth be referred to as "T1".

Construction of a 3D Modeled Surface at the Center of the Seam
As stated above, once the borehole intervals that intercept the mineralized body have been defined, a surface at the center of this body has to be created, from where the interpolation point mesh will be created.As noted above, the surface does not need to be located exactly, as the software will adjust its final position according to the real values of the borehole intercepts and the body.
The easiest way to proceed is to work section by section (Figure 4). Figure 4A shows a cross section of boreholes depicted in white whereas the intercepts with the body are depicted in red.The green line connects the centers of each red line section.Figure 4B simultaneously shows several sections and their corresponding centerlines, which are triangulated to generate a surface.Figure 4C shows the final calculated surface, which will henceforth be referred to as "T1".

Generation of an Equally Spaced Cloud of Points in the Seam Centered Surface
In order to have equally spaced interpolation points over the T1 surface, T1 will be remeshed to generate T2, with a point density high enough to reflect geology variation.A high point density will not alter the resource calculations, as interpolation points too close to each other will display the same results as if they were close to the intercepts.This cloud of points will be referred to as NPS from now onwards.The software will automatically cut parallel sections of T1 at a defined distance (5 m in this example, as the borehole density is 25 × 25 and 5 intermediate points between boreholes have been considered) in a selected direction (the X axis in this example).Inside each line generated by each cut, points are defined by a certain distance, user-selectable (5 m again in this example) on the perpendicular axis (Y in this case).Note that 5 m is a suitable value in the case of the "ZONA M" ore body.This value has to be adjusted according to each ore body size and morphology so as to obtain the necessary accuracy for mining the resource, see [7].
Once the evenly spaced cloud of points is generated (Cartesian mesh), the software will create new points at the cross of the diagonals connecting any four points of each generated square (Figure 5C). Figure 5B shows the outcome of this operation.Finally, the algorithm will create triangles that mimic the seam-centered surface.This triangulated surface will be called T2.
All the information is stored in a database (BDS), previously created in the database.

Generation of an Equally Spaced Cloud of Points in the Seam Centered Surface
In order to have equally spaced interpolation points over the T1 surface, T1 will be remeshed to generate T2, with a point density high enough to reflect geology variation.A high point density will not alter the resource calculations, as interpolation points too close to each other will display the same results as if they were close to the intercepts.This cloud of points will be referred to as NPS from now onwards.The software will automatically cut parallel sections of T1 at a defined distance (5 m in this example, as the borehole density is 25 × 25 and 5 intermediate points between boreholes have been considered) in a selected direction (the X axis in this example).Inside each line generated by each cut, points are defined by a certain distance, user-selectable (5 m again in this example) on the perpendicular axis (Y in this case).Note that 5 m is a suitable value in the case of the "ZONA M" ore body.This value has to be adjusted according to each ore body size and morphology so as to obtain the necessary accuracy for mining the resource, see [7].
Once the evenly spaced cloud of points is generated (Cartesian mesh), the software will create new points at the cross of the diagonals connecting any four points of each generated square (Figure 5C). Figure 5B shows the outcome of this operation.Finally, the algorithm will create triangles that mimic the seam-centered surface.This triangulated surface will be called T2.

Generation of an Equally Spaced Cloud of Points in the Seam Centered Surface
In order to have equally spaced interpolation points over the T1 surface, T1 will be remeshed to generate T2, with a point density high enough to reflect geology variation.A high point density will not alter the resource calculations, as interpolation points too close to each other will display the same results as if they were close to the intercepts.This cloud of points will be referred to as NPS from now onwards.The software will automatically cut parallel sections of T1 at a defined distance (5 m in this example, as the borehole density is 25 × 25 and 5 intermediate points between boreholes have been considered) in a selected direction (the X axis in this example).Inside each line generated by each cut, points are defined by a certain distance, user-selectable (5 m again in this example) on the perpendicular axis (Y in this case).Note that 5 m is a suitable value in the case of the "ZONA M" ore body.This value has to be adjusted according to each ore body size and morphology so as to obtain the necessary accuracy for mining the resource, see [7].
Once the evenly spaced cloud of points is generated (Cartesian mesh), the software will create new points at the cross of the diagonals connecting any four points of each generated square (Figure 5C). Figure 5B shows the outcome of this operation.Finally, the algorithm will create triangles that mimic the seam-centered surface.This triangulated surface will be called T2.
All the information is stored in a database (BDS), previously created in the database.All the information is stored in a database (BDS), previously created in the database.

Calculation of the Intercept Values of the Surveys Crossing the Seam
Before the interpolation of grades and thicknesses at each point of the NPS, the economical interval in each borehole (intercept) must be calculated.
Cutoff grade.This is defined as the minimum grade that allows economical exploitation of the ore body, see [19].In the example that is being used, the "Zona M" ore body of the Carlés deposit hosts three recoverable elements: Au, Cu, and Ag.A new element called "Au_e" element or "Au grade equivalent" can be defined, where an accumulated value of the three elements is stored according to a price-weighting index scheme.In order to simplify the calculation, an average metal recovery value was assumed, understood as the final percentage of each metal that gives a monetary income, deducting the concentration plant recovery, selling costs, concentration transport, penalties, besides other factors.
Minimum thickness.A minimum value of seam thickness must be considered in order to reflect the economic feasibility of the exploitation.It depends on the rock type, the seam dip and the mining method (room and pillar, sublevel stopping, cut and fill).A thickness of 2 m has been considered in this example, according to the mining machinery size that can be used in this deposit mining.
Mining overbreak.A lateral dilution has been defined due to unwanted results of blasting and hole deviations.From our experience with this deposit, 50 cm of overbreak (overexcavation, 25 cm each side) is defined, a value that can be user-selectable Maximum thickness of waste rock.A usual situation that is found when the intercept of the borehole and the ore body is calculated is the appearance of inclusions, intervals of waste rock inside the ore body and, more commonly, mineralization below the cutoff grade value.The treatment for this type of situation is achieved by defining a maximum waste thickness value whereby, if the waste thickness value exceeds the maximum, it will not be mined (it will remain as a pillar).If the waste thickness value is calculated to be below the maximum and the average grade is higher than the cutoff grade, it will be considered as ore; however, if the average grade is below the cutoff grade, it will be considered as waste rock and will not be included in the intercept.This maximum waste rock thickness will be referred to as the real thickness of the seam, measured perpendicular to the seam.
The intercept calculation procedure is described below, both graphically and analytically.
In order to explain the procedure of the three intersection types, as well as the interpolation over each point of the NPS, we will refer to four theoretical boreholes: S1, S2, S3, and S4, as represented in Figure 6.The algorithm calculates three different intercepts (denoted as A, B, and C), for each borehole crossing the seam, as follows: A: Geological Intercept.This includes only those samples with grades higher than the cutoff grade.If there are no samples higher than that value, the highest-grade sample will be defined as the geological intercept.For instance, borehole S2 shows six samples 1 m in size, but only five will be used when considering a cutoff grade of 2 g Au/t, as can be seen in Figure 4 as per the coding colors.
For those samples with grades lower than the cutoff grade being inserted between samples that are above the cutoff grade, the "maximum thickness of waste rock" criteria must be taken into account, as referred to above.
B: Minimum Thickness Intercept.This will be defined as the geological intercept widened as necessary to complete a previously defined minimum real thickness.If the real thickness of the geological intercept is equal to or greater than the minimum thickness value, both intercepts will be the same.
By looking at borehole S2 in Figure 6, and considering 2 m as the minimum thickness value, it can be seen that the geological intercept already has a real thickness of over 2 m, so the minimum thickness intercept (B) will be equal to the geological intercept (A).
If the intercept is less than the minimum thickness, the algorithm will add in samples from the higher grade side.Once the limits are checked, it will re-calculate whether new samples need to be added of a higher cutoff on each side of the intercept.C: Mining Overbreak Intercept.This will be defined as the minimum real thickness intercept widened at both sides by a distance defined as over excavation.In this example 50 cm will be used (25 cm on each side).This is an important parameter to consider, as it will simulate the mineral dilution that always occurs in common mining procedures that involve excavation.
this type of situation is achieved by defining a maximum waste thickness value whereby, if the waste thickness value exceeds the maximum, it will not be mined (it will remain as a pillar).If the waste thickness value is calculated to be below the maximum and the average grade is higher than the cutoff grade, it will be considered as ore; however, if the average grade is below the cutoff grade, it will be considered as waste rock and will not be included in the intercept.This maximum waste rock thickness will be referred to as the real thickness of the seam, measured perpendicular to the seam.
The intercept calculation procedure is described below, both graphically and analytically.
In order to explain the procedure of the three intersection types, as well as the interpolation over each point of the NPS, we will refer to four theoretical boreholes: S1, S2, S3, and S4, as represented in Figure 6.The algorithm calculates three different intercepts (denoted as A, B, and C), for each borehole crossing the seam, as follows:  All three intercepts are calculated referring to real or true thickness; therefore the angle at which the borehole and the seam intercept each other must be known.If the borehole is perpendicular to the seam, this angle is 90 • and real and apparent thickness will be equal.This angle, and the intercept point coordinates, is calculated automatically by the algorithm.The intercept point will be the point where the line that represents the borehole cuts one of the triangles of the T2 surface generated from the NPS.The angle will be the one established by the plane of the intersected triangle and the segment representing the borehole.
Nevertheless, the experience of the technician is extremely important in order to define which samples to include in the body.The software allows samples to be accepted or discarded manually in order to refine the calculations.Angles defined by the intersection of boreholes and the seam shaped deposit are calculated at each intersection and can be manually adjusted according to geological knowledge if need be.In this way, possible inaccuracies in apparent thickness calculations are limited.
Interpolation parameters.In this study the Inverse Distance Weighted (IDW) interpolation method will be used; other methods could have been considered: nearest neighbor or kriging, among others.In order to use IDW, a maximum influence distance must be used, calculated through a geostatistical study, independent of the interpolation method shown in this paper and detailed below, as well as a power value.In this case the distance is fixed at 50 m and the power is 3 (Figure 7 shows the input window for these parameters).These values can be different for each intercept, if necessary.
Once the calculation parameters are defined, the software executes the algorithm, which resolves the three intercept types in each borehole that fulfills the defined conditions.This kind of deposit is often high-grade centered, with low-grade values toward the deposit limits.When several high-grade intervals exist over the cutoff grade, the limits for each and every one of the three types of intercepts can be adjusted manually.The calculation results (Figure 9) are stored in a database (BDS) that includes three records for each borehole, corresponding to each interpolation type.The calculation results (Figure 9) are stored in a database (BDS) that includes three records for each borehole, corresponding to each interpolation type.

Geostatistical Study of the Intercept Values
As with block model estimation, an influence distance or neighborhood search must be defined.This is necessary in order to consider the quantity and proximity of intercepts to be used in the interpolation and to define the parameters of the classification of the resources or reserves.This is achieved by using a variographic study, which can be performed using a variety of software to which the data must be exported.Several codes are available, both freeware and commercial.
In this example, the data were exported to Stanford Geostatistical Modeling Software (SGEMS), [20], where the variogram was modeled and the main direction and influence distance was found.RecMin can easily export both the composite and/or the intercept data to make this calculation, in both Gslib and Ascii format.
The main direction of the mineralization, as per variogram results, is indicated in graphics by a line parallel to the direction that fits between two of the vertices of the T2 surface (Figure 10).This line is saved as a file and then later selected in the software interpolation assistant in order to be considered as a continuity guide during calculation.
The anisotropy of the variogram is reflected in the different radii of the search ellipsoid.In this case, it is a planar ellipse, as it follows a seam.Both radii of the ellipse, the long and the short, are introduced considering a ratio between them.In this example the ratio is 1:2, with 2 being the secondary direction.Figure 10 shows the shape of the ellipse for a unique theoretical borehole, crossing the seam at the center of the colored sections, as well as the main direction of the mineralization (green line), which follows the main radius of the ellipse.

Geostatistical Study of the Intercept Values
As with block model estimation, an influence distance or neighborhood search must be defined.This is necessary in order to consider the quantity and proximity of intercepts to be used in the interpolation and to define the parameters of the classification of the resources or reserves.This is achieved by using a variographic study, which can be performed using a variety of software to which the data must be exported.Several codes are available, both freeware and commercial.
In this example, the data were exported to Stanford Geostatistical Modeling Software (SGEMS), [20], where the variogram was modeled and the main direction and influence distance was found.RecMin can easily export both the composite and/or the intercept data to make this calculation, in both Gslib and Ascii format.
The main direction of the mineralization, as per variogram results, is indicated in graphics by a line parallel to the direction that fits between two of the vertices of the T2 surface (Figure 10).This line is saved as a file and then later selected in the software interpolation assistant in order to be considered as a continuity guide during calculation.
The anisotropy of the variogram is reflected in the different radii of the search ellipsoid.In this case, it is a planar ellipse, as it follows a seam.Both radii of the ellipse, the long and the short, are introduced considering a ratio between them.In this example the ratio is 1:2, with 2 being the secondary direction.Figure 10 shows the shape of the ellipse for a unique theoretical borehole, crossing the seam at the center of the colored sections, as well as the main direction of the mineralization (green line), which follows the main radius of the ellipse.If we connect the ends of the three segments on each corresponding side, two new triangles appear; these triangles, together with the three segments, define a volume that is called the "calculation unit" (Figure 11a).If we connect the ends of the three segments on each corresponding side, two new triangles appear; these triangles, together with the three segments, define a volume that is called the "calculation unit" (Figure 11a).The volume is calculated by dividing each calculation unit (irregular octahedron) in three irregular tetrahedra, as can be seen in Figure 9b.As the coordinates of each of the four vertices of the irregular tetrahedron are known, its volume can be calculated as shown in Equation (10), where xi, yi, zi are the vertex coordinates: It must be taken into account that each pair of triangles that share an edge will also share the two segments generated at each vertex edge and the corresponding faces of those two triangles will also be the same.In this way, all calculation units match and fill the entire volume that represents the seam (Figure 12).The volume is calculated by dividing each calculation unit (irregular octahedron) in three irregular tetrahedra, as can be seen in Figure 9b.As the coordinates of each of the four vertices of the irregular tetrahedron are known, its volume can be calculated as shown in Equation (10), where x i , y i , z i are the vertex coordinates: It must be taken into account that each pair of triangles that share an edge will also share the two segments generated at each vertex edge and the corresponding faces of those two triangles will also be the same.In this way, all calculation units match and fill the entire volume that represents the seam (Figure 12).Finally, the weight will be obtained by multiplying the calculated volume by the seam density for each one of the three calculation types.
Before initiating the calculation, all the parameters are introduced into the software (Figure 13).Finally, the weight will be obtained by multiplying the calculated volume by the seam density for each one of the three calculation types.
Before initiating the calculation, all the parameters are introduced into the software (Figure 13).Finally, the weight will be obtained by multiplying the calculated volume by the seam density for each one of the three calculation types.
Before initiating the calculation, all the parameters are introduced into the software (Figure 13).

Interpolation of the Positions of the Seam Surface Center
As we have seen in the calculation of intersections, once the calculations are complete and depending on the cutoff grade, minimum thickness or overbreak dilution, the intercept center does not have to match the theoretical center where the borehole cuts the seam center (T1).Differences at these locations will not affect a volume and grade calculation.Nevertheless, it will be of significance if the calculations are going to be used in the planning or exploitation mining phases; the way to correct these positions is by using position interpolation.In this way, the software relocates the center of the seam, moving it depending on how the nearest intercept centers have been moved, and the information is then stored in database fields named MOVX, MOVY, and MOVZ.This procedure relocates the NPS interpolation points by applying the inverse distance weighting method with a certain user-selectable power.
The software can be instructed to store the new point mesh positions in a new table so they can be reused in a subsequent calculation and to obtain graphical representations (Figure 14).

Interpolation of the Positions of the Seam Surface Center
As we have seen in the calculation of intersections, once the calculations are complete and depending on the cutoff grade, minimum thickness or overbreak dilution, the intercept center does not have to match the theoretical center where the borehole cuts the seam center (T1).Differences at these locations will not affect a volume and grade calculation.Nevertheless, it will be of significance if the calculations are going to be used in the planning or exploitation mining phases; the way to correct these positions is by using position interpolation.In this way, the software relocates the center of the seam, moving it depending on how the nearest intercept centers have been moved, and the information is then stored in database fields named MOVX, MOVY, and MOVZ.This procedure relocates the NPS interpolation points by applying the inverse distance weighting method with a certain user-selectable power.
The software can be instructed to store the new point mesh positions in a new table so they can be reused in a subsequent calculation and to obtain graphical representations (Figure 14).
The change in position of the seam center will be the same in the case of "Minimum thickness" and "Mining" interpolations, as they have the same center; in contrast, it will be different in the case of the "Geological" type, so in this case it does not need to be the same, as explained previously.

Definition of the Estimation Categories and Generation of the Calculation Units
As in any other resource estimation method, the interpolation level of confidence must be defined.This is done through the use of the so-called categories.The algorithm uses up to four categories (measured, indicated, inferred or proven, and probable).The position of samples in the octant can be used if necessary, avoiding and/or adjusting calculations according to the number of samples present in opposite octants surrounding the calculation point.In the example shown here, four numerical categories are used: • Category 1: calculation units with an intercept at less than 10 m The change in position of the seam center will be the same in the case of "Minimum thickness" and "Mining" interpolations, as they have the same center; in contrast, it will be different in the case of the "Geological" type, so in this case it does not need to be the same, as explained previously.

Definition of the Estimation Categories and Generation of the Calculation Units
As in any other resource estimation method, the interpolation level of confidence must be defined.This is done through the use of the so-called categories.The algorithm uses up to four categories (measured, indicated, inferred or proven, and probable).The position of samples in the octant can be used if necessary, avoiding and/or adjusting calculations according to the number of samples present in opposite octants surrounding the calculation point.In the example shown here, four numerical categories are used:

•
Category 1: calculation units with an intercept at less than 10 m • Category 2: calculation units with two intercepts at less than 20 m • Category 3: calculation units with two intercepts at less than 30 m • Category 4: calculation units with one intercept at less than 40 m In these category calculations, the distance reduction defined by the main and secondary axis of the ellipsoid is also taken into account.

Calculation Units Generation and Categories Definition
Some parts of the calculation method are based on the generation of provisional data sheets created through the addition of data queries.This method allows fast operations across the entire data set.
The SQL language was also used to make the selections and differentiate, quickly and with low coding between the different calculation types.
The generated information was stored in a relational database (BDT2), from where it is retrieved by the drawing module.Figure 15 shows the last tab of the calculation window, where the total amounts corresponding to a certain cutoff grade can be seen, as well as the total for each category.In these category calculations, the distance reduction defined by the main and secondary axis of the ellipsoid is also taken into account.

Calculation Units Generation and Categories Definition
Some parts of the calculation method are based on the generation of provisional data sheets created through the addition of data queries.This method allows fast operations across the entire data set.
The SQL language was also used to make the selections and differentiate, quickly and with low coding effort, between the different calculation types.
The generated information was stored in a relational database (BDT2), from where it is retrieved by the drawing module.Figure 15 shows the last tab of the calculation window, where the total amounts corresponding to a certain cutoff grade can be seen, as well as the total for each category.The influence by percentage of each borehole of the total amount of weight and interpolated elements can also be obtained (Table 3).This allows for the checking of the results, specifically searching for balanced results and avoiding overweighed boreholes.The main results can be guided by a very low number of boreholes, which can result in a significant variation if sample mistakes are present.Figure 16 shows SQL guided results analysis, which allows for powerful, quick, and easy selection of data.Figure 17A shows the 3D results of the SQL selection and the associated legend colors.Please notice that the seam width variation is adequately simulated, as shown in Figure 17B.The influence by percentage of each borehole of the total amount of weight and interpolated elements can also be obtained (Table 3).This allows for the checking of the results, specifically searching for balanced results and avoiding overweighed boreholes.The main results can be guided by a very low number of boreholes, which can result in a significant variation if sample mistakes are present.Figure 16 shows SQL guided results analysis, which allows for powerful, quick, and easy selection of data.Figure 17A shows the 3D results of the SQL selection and the associated legend colors.Please notice that the seam width variation is adequately simulated, as shown in Figure 17B.

Method Suitability
This method has been called the "Pentahedral Method" due to the fact that the calculation unit has a shape close to a pentahedron: two triangular faces connected by three four-sided faces.Nevertheless, it is clear that it is not a pentahedron, but rather an octahedron, as the faces with four sides could have their vertices on different planes; that is, each of them is made up of two triangular faces (Figure 11A).
This method is a calculation alternative in the case of these seam-shaped mineral structures, as it can work in three dimensions and be completely computer-programmable.
It is suited to seams with extremely pronounced folds, as seen in Figure 18.

Method Suitability
This method has been called the "Pentahedral Method" due to the fact that the calculation unit has a shape close to a pentahedron: two triangular faces connected by three four-sided faces.Nevertheless, it is clear that it is not a pentahedron, but rather an octahedron, as the faces with four sides could have their vertices on different planes; that is, each of them is made up of two triangular faces (Figure 11A).
This method is a calculation alternative in the case of these seam-shaped mineral structures, as it can work in three dimensions and be completely computer-programmable.
It is suited to seams with extremely pronounced folds, as seen in Figure 18.

Method Suitability
This method has been called the "Pentahedral Method" due to the fact that the calculation unit has a shape close to a pentahedron: two triangular faces connected by three four-sided faces.Nevertheless, it is clear that it is not a pentahedron, but rather an octahedron, as the faces with four sides could have their vertices on different planes; that is, each of them is made up of two triangular faces (Figure 11A).
This method is a calculation alternative in the case of these seam-shaped mineral structures, as it can work in three dimensions and be completely computer-programmable.
It is suited to seams with extremely pronounced folds, as seen in Figure 18.Compared to classic block models the pentahedral method offers several advantages, such as:  Any change in the calculation parameters does not need the calculation units to be redefined.

•
When defining the calculation geometry units, these can be separately prepared for subsequent exploitation and scheduling, even exporting them to other software.

•
It provides a fairly accurate representation of the thickness of the seam-shaped structure even in the case of low thickness seams, which is extremely difficult to achieve with traditional block models (Figure 19).

•
It interpolates intersections, not samples taken from boreholes.

•
Geometrical results closely follow the deposit limits, as opposed to block models that create a stair-shaped profile (see Figure 19).

•
It includes the possible internal dilution, as well the lateral dilution due to mining overbreak.

•
A minimum thickness can be defined according to mining technical conditions.
Minerals 2017, 7, 72 19 of 21 Compared to classic block models the pentahedral method offers several advantages, such as: • Any change in the calculation parameters does not need the calculation units to be redefined.• When defining the calculation geometry units, these can be separately prepared for subsequent exploitation and scheduling, even exporting them to other software.

•
It provides a fairly accurate representation of the thickness of the seam-shaped structure even in the case of low thickness seams, which is extremely difficult to achieve with traditional block models.(Figure 19).

•
It interpolates intersections, not samples taken from boreholes.• Geometrical results closely follow the deposit limits, as opposed to block models that create a stair-shaped profile (see Figure 19).• It includes the possible internal dilution, as well as the lateral dilution due to mining overbreak.• A minimum thickness can be defined according to mining technical conditions.
This method could be also used to model larger ore deposits or even regional scale geological seams, such as the ones shown in [21].

Discussion and Conclusions
This paper shows that, in the case of mineral ore bodies that are seam-shaped, using the new "Pentahedral Method" it is possible to obtain an evaluation of the reserves and resources with a higher degree of accuracy than with the classic block method.
Compared to other methods such as the polygon method, the triangle method, or the section one, this methodology is more advanced, allowing interpolation with the intercepts.This method can be applied to all types of seams, even those with extremely pronounced folds, and it allows work to be conducted in 3D and to be a good representation of the ore body to be obtained.Figure 20 is a U3D model that can be embedded in the pdf of the paper, which shows interactively what has been explained throughout this paper.This method could be also used to model larger ore deposits or even regional scale geological seams, such as the ones shown in [21].

Discussion and Conclusions
This paper shows that, in the case of mineral ore bodies that are seam-shaped, using the new "Pentahedral Method" it is possible to obtain an evaluation of the reserves and resources with a higher degree of accuracy than with the classic block method.
Compared to other methods such as the polygon method, the triangle method, or the section one, this methodology is more advanced, allowing interpolation with the intercepts.This method can be applied to all types of seams, even those with extremely pronounced folds, and it allows work to be conducted in 3D and to be a good representation of the ore body to be obtained.Figure 20 is a U3D model that can be embedded in the pdf of the paper, which shows interactively what has been explained throughout this paper.
This method also offers several additional advantages compared to the classic block model concept, mainly in its three-dimensional representation capabilities.Regarding interpolation, it also allows any method to be used, even kriging, as well as the use of search ellipsoids.This method also offers several additional advantages compared to the classic block model concept, mainly in its three-dimensional representation capabilities.Regarding interpolation, it also allows any method to be used, even kriging, as well as the use of search ellipsoids.
The use of SQL databases to store the results allows reports and results analysis to be issued very quickly.
The introduction of calculation variables such as "overbreak," "dilution," or "minimum thickness" gives more information from the resources data to the technician in charge of the reserves definition.These variables are related to the mining process and will simplify the reserve calculation as it will eliminate some of the resources that are not economically exploitable.
These results can be used for mine planning and design where final reserves will be defined, and have already been included in some versions of the Recmin software.At this stage the method has not been prepared to be used with other software, although it could be adapted to generate data and share it through simple ASCII files where block data are coded.
The method has passed the preliminary exam of the international patent PCT/ES03/00117, achieving positive results; the research is recognized as innovative, it implies research activity and it can be applied in the industry.

Figure 2 .
Figure 2. (A) View of the stratigraphic control of the skarn development at the Carlés skarn.The monzogranite intrusion is also controlled by the sedimentary layers.(B) Close-up view of stratigraphically controlled Carlés mineralized banded skarn.Bands are dominated by pyroxenemagnetite alternate with garnet-amphibole-rich bands.

Figure 2 .
Figure 2. (A) View of the stratigraphic control of the skarn development at the Carlés skarn.The monzogranite intrusion is also controlled by the sedimentary layers; (B) Close-up view of stratigraphically controlled Carlés mineralized banded skarn.Bands are dominated by pyroxene-magnetite alternate with garnet-amphibole-rich bands.

Figure 3 .
Figure 3. Example of a mineralized area in an intercepting borehole.

Figure 3 .
Figure 3. Example of a mineralized area in an intercepting borehole.

Figure 4 .
Figure 4. T1 surface generation.(A) A cross section of boreholes is depicted in white and the intercepts with the body are depicted in red.The green line connects the centers of each red line section.(B) Sections and their corresponding centerlines, which are triangulated to generate a surface.(C) Calculated surface ("T1").

Figure 5 .
Figure 5. Example of base surface T1 (left (A)) converted to calculation Surface T2 (right (B)).Detail of meshing points in sections and in crossing diagonals of any generated square (C).

Figure 4 .
Figure 4. T1 surface generation.(A) A cross section of boreholes is depicted in white and the intercepts with the body are depicted in red.The green line connects the centers of each red line section; (B) Sections and their corresponding centerlines, which are triangulated to generate a surface.(C) Calculated surface ("T1").

Figure 4 .
Figure 4. T1 surface generation.(A) A cross section of boreholes is depicted in white and the intercepts with the body are depicted in red.The green line connects the centers of each red line section.(B) Sections and their corresponding centerlines, which are triangulated to generate a surface.(C) Calculated surface ("T1").

Figure 5 .
Figure 5. Example of base surface T1 (left (A)) converted to calculation Surface T2 (right (B)).Detail of meshing points in sections and in crossing diagonals of any generated square (C).

Figure 5 .
Figure 5. Example of base surface T1 (left (A)) converted to calculation Surface T2 (right (B)).Detail of meshing points in sections and in crossing diagonals of any generated square (C).

Figure 6 .
Figure 6.Intercept calculation of four boreholes (S1 to S4) through a seam-shaped mineral ore body.A (geology intercept), B (minimum thickness intercept), C (mining overbreak intercept); U, V, and W are points of the NPS; i1-i4 are the lithology intercept interval; d1-d4 distances between the intercepts center and W point.

Figure 6 .
Figure 6.Intercept calculation of four boreholes (S1 to S4) through a seam-shaped mineral ore body.A (geology intercept), B (minimum thickness intercept), C (mining overbreak intercept); U, V, and W are points of the NPS; i1-i4 are the lithology intercept interval; d1-d4 distances between the intercepts center and W point.

Figure 8
Figure 8 is an example of the calculation window.Intercepts are shown as graph bars: red (geological intercept), yellow (minimum thickness intercept) and white (mining overbreak intercept).The x-axis shows samples ("from" and "to" values) and the y-axis shows the different grade values in ppb (parts per billion).

Figure 8
Figure 8 is an example of the calculation window.Intercepts are shown as graph bars: red (geological intercept), yellow (minimum thickness intercept) and white (mining overbreak intercept).The x-axis shows samples ("from" and "to" values) and the y-axis shows the different grade values in ppb (parts per billion).

Figure 8
Figure 8 is an example of the calculation window.Intercepts are shown as graph bars: red (geological intercept), yellow (minimum thickness intercept) and white (mining overbreak intercept).The x-axis shows samples ("from" and "to" values) and the y-axis shows the different grade values in ppb (parts per billion).

Figure 9 .
Figure 9. Table of calculated intercepts calculated in RecMin.Results for several boreholes.First column is the name of the borehole (down hole).Each hole has three data rows showing the intercept data.

Figure 9 .
Figure 9. Table of calculated intercepts calculated in RecMin.Results for several boreholes.First column is the name of the borehole (down hole).Each hole has three data rows showing the intercept data.

Figure 11 .
Figure 11.Calculation unit as used in the pentahedral method (A).Volume calculation (B) and calculation graphic exploded (C).

Figure 11 .
Figure 11.Calculation unit as used in the pentahedral method (A).Volume calculation (B) and calculation graphic exploded (C).

Figure 12 .
Figure 12.Mineral seam filled with calculation units, as automatically obtained by the software.Boreholes are colored by lithology.The legend refers to grade over seam (ppb).

Figure 13 .
Figure 13.Calculation parameters in the interpolation of the pentahedral method.

Figure 12 .
Figure 12.Mineral seam filled with calculation units, as automatically obtained by the software.Boreholes are colored by lithology.The legend refers to grade over seam (ppb).

Minerals 2017, 7 , 72 14 of 21 Figure 12 .
Figure 12.Mineral seam filled with calculation units, as automatically obtained by the software.Boreholes are colored by lithology.The legend refers to grade over seam (ppb).

Figure 13 .
Figure 13.Calculation parameters in the interpolation of the pentahedral method.

Figure 13 .
Figure 13.Calculation parameters in the interpolation of the pentahedral method.

Figure 14 .
Figure 14.Example of economical seam center (blue line) vs. theoretical seam center (green line).

Figure 14 .
Figure 14.Example of economical seam center (blue line) vs. theoretical seam center (green line).

Figure 15 .
Figure 15.Numerical results of the "Pentahedral Method," classified by intercept type and category.

Figure 15 .
Figure 15.Numerical results of the "Pentahedral Method," classified by intercept type and category.

Figure 16 .
Figure 16.Assistant module used to open the pentahedral files and obtain a 3D visualization, as per SQL database conditions restrictions.Figure 16.Assistant module used to open the pentahedral files and obtain a 3D visualization, as per SQL database conditions restrictions.

Figure 16 .
Figure 16.Assistant module used to open the pentahedral files and obtain a 3D visualization, as per SQL database conditions restrictions.Figure 16.Assistant module used to open the pentahedral files and obtain a 3D visualization, as per SQL database conditions restrictions.

Figure 17 .
Figure 17.3D view, pentahedral filled seam representing orebody "Zona M" of the Carlés deposit (A).Detail of the variable width of the calculated seam (B).The legend shows the Au grade, expressed in ppb.

Figure 17 .
Figure 17.3D view, pentahedral filled seam representing orebody "Zona M" of the Carlés deposit (A).Detail of the variable width of the calculated seam (B).The legend shows the Au grade, expressed in ppb.

Minerals 2017, 7 , 72 18 of 21 Figure 17 .
Figure 17.3D view, pentahedral filled seam representing orebody "Zona M" of the Carlés deposit (A).Detail of the variable width of the calculated seam (B).The legend shows the Au grade, expressed in ppb.

Figure 19 .
Figure 19.Graphic comparison between the block model (conventional method) and the "Pentahedral Method".Left (A): plane section of the blocks that define the seam colored according to grade and the same calculated seam.Right (B): sections of the three different types of calculation done in the same seam.

Figure 19 .
Figure 19.Graphic comparison between the block model (conventional method) and the "Pentahedral Method".Left (A): plane section of the blocks that define the seam colored according to grade and the same calculated seam; Right (B): sections of the three different types of calculation done in the same seam.

Figure 20 .
Figure 20.3D model of the calculated seam colored by grade, including in the user selectable layers the topography and boreholes colored according to geology and sample grading.

Table 2 .
Calculation breakdown and interpolation of an NPS triangle.

Table 2 .
Calculation breakdown and interpolation of an NPS triangle.

Table 3 .
Influence percentage by borehole.

Table 3 .
Influence percentage by borehole.