APPLICATION OF AN INVERSE DISTANCE WEIGHTED ANISOTROPIC METHOD FOR ROCK QUALITY DESIGNATION DISTRIBUTION IN EASTERN KAHANG DEPOSIT, CENTRAL IRAN

Estimation of Rock Quality Designation (RQD) distribution with low value of error is crucial for mining excavation based on geomechanical data. This paper proposes an application of an Inverse Distance Weighted Anisotropic Method (IDWAM) based on combined variograms to construct RQD block model in Kahang Copper-Molybdenum (Cu-Mo) porphyry deposit, Central Iran, using subsurface data. To do this, an appropriate voxel size was calculated and then RQD variography was carried in horizontal and vertical directions. In addition, Deere and Miller rock classification was used to classify RQD block model for a final open pit interpretation. Finally, correlation between results for RQD estimated derived via the IDWAM and raw data was carried out using jackknife method. The results obtained by the combination of IDWAM and experimental variogram showed that the excellent RQD zones are located in the central and North Western parts of this area.


Introduction
Assessment of Rock Quality Designation (RQD) is a crucial aspect for mineral excavation, resource modelling and mine planning with huge cost implications for the design and the mining of each block such as ore or waste tonnage [1][2][3]. These are calculated using the dimension and density of each block. The results are assessed to identify a final pit slope angle and consequent pit stripping ratio [4][5][6][7]. Numerical models in geosciences have been created and consequently utilised to better interpret the variability of geological parameters such as lithology, ore-type, alteration and mineralogy or for a better understanding of the different attributes such as density, rock mass characterisation and RQD [8][9][10]. However, the classical statistical methods for delineation of populations from a background level would be for example, a histogram analysis, box plot, summation of mean and standard deviation coefficients and median. These are not considered overly accurate due to the fact that these statistical methods consider only the frequency distribution of information while not paying attention to the spatial variability [11,12]. In other words, the classical statistical plots (i.e., histograms) are based on the data abundance distribution and cannot quantify the spatial positions of parameters such as RQD [13][14][15][16]. As a result, numerical modelling of rock characteristics is a difficult task and requires 3D modelling for better interpretation of the problems found in a mining operation such as rock discontinuities, planar failure, circular failure, wedge failure and toppling failure [17,18].
The earliest model regarding the quantitative description of in-situ block size distribution (IBSD) was the Rock Quality Designation, which was proposed by [19]. Priest and Hudson (1976) applied the RQD method to scanline survey data with respect to an analytical relation between RQD and the discontinuity frequency resulted from a scanline survey [20]. A borehole or a scanline are by their nature one dimensional. As a result, RQD values calculated using bore hole data or a scanline survey are influenced by the orientation in which the measurements are taken (horizontal or vertical) so the method does not consider calculation for the other direction [21,22]. In order to overcome the dependence of RQD on orientation, Kazi and Sen (1985) proposed the use of the Volumetric Rock Quality Designation (V. RQD) which is a three-dimensional parameter [23]. This encloses the proportion of the volume of intact matrix rock blocks, equal to or higher than 0.001 m3 in size, which can be associated with the average volume of a matrix block and the number of matrix blocks per m 3 . To do this, the V.RQD is calculated by summation of the volumes of intact blocks divided by the total rock mass volume which is expressed as a percentage. However, the proposed model is limited to the estimation of the average block volume rather than the IBSD [19,23]. Palmstrom (1985) proposed different empirical equations to link Volumetric Discontinuity Count for RQD (known as Jv RQD) data and linear fracture frequency [24]. He suggested that there is a correlation between the in-situ block size and Jv is represented in a figure incorporating various measurements of the block size or degree of jointing (e.g., density of joints, RQD, block volume and joint spacing). However, this model can only estimate a rough upper and lower range of block sizes and therefore has restricted practical applications [25]. Şen and Eissa (1992) derived values for Jv for RQD and block volumes of different shapes such as bars, plates and or prisms quantity [23]. This model which was proposed by Şen and Eissa (1992) provides a simple tool for rock engineers without the need for recourse to theoretical calculations. However, the block volume in this model is given in terms of average block size so it cannot describe the block size distribution [26].
Techniques for selecting a theoretical function to describe rock mass characteristics based on RQD data are unsatisfactory, as mentioned above. As a result, an approach to introduce an analytical model to help achieve 2D and/or 3D maps for interpreting the distributions of measured RQD is needed. To do this, geostatistical methods can be used for interpolation and estimation of different regional variables (RQD in this scenario). Employment of a proper estimation method regarding geometry and geological properties as well as drilling patterns for different ore deposits is a fundamental issue in resource estimation [27][28][29][30][31]. Linear/ non-linear Kriging methods and Inverse Distance Weighted (IDW) have been widely utilised for block modelling in mineral exploration [32][33][34]. This paper aims to introduce an application of ''Inverse Distance Weighted Anisotropic Method (IDWAM)'' to create an RQD block model in Kahang Cu-Mo porphyry deposit, central Iran. The obtained RQD block model is validated via the Deere and Miller rock classification (Table 1) to suggest an optimised final pit slope.

Geological setting of the Kahang Cu-Mo porphyry deposit
The Kahang deposit is located about 73 km NE of Isfahan in central Iran. The deposit is situated in the Cenozoic Urumieh-Dokhtar magmatic belt extending from NW to SE Iran, as depicted in Fig.  1 [35][36][37]. This deposit is mainly composed of Eocene volcanic-pyroclastic rocks, which were intruded by quartzmonzonite, monzogranite to dioritic intrusions in Oligo-Miocene rocks (Fig. 1). The extrusive rocks include tuffs, breccias and lavas. The existing lavas in this area have dacitic to andesitic composition. The main criteria aspects for determining detection of mineralised zones are their index ores minerals. Index ores which consist of chalcopyrite and pyrite for hypogene zone, . Furthermore, chalcocite, bornite and covellite for are index minerals for supergene enrichment zone. and malachiteMalachite, azurite, tenorite and cuprite occur for within oxidation zone [38]. Studies Investigation of the pattern of zonation in the eastern part of the Kahang deposit show indicates that the most significant major mineralisation (in terms of ore zone size and geometry) is hypogene zone which containing includes a high percentage volumes of chalcopyrite accompanied byand pyrite. The major main alteration zones of are potassic, phyllic, argillic and propylitic types were accompanied by thewith vein to veinlets fillings of quartz, quartz-magnetite and Fe-hydroxides [39][40].

Methodology
In this deposit, 14979 core samples were determined from 42 boreholes with total depth of ~22,000 m in the deposit for RQD analysis. Secondly, the subsurface data, including coordinates of drillcores, azimuth and dip (orientation) and RQD values (measured) was generated via RockWorksTM software package. Then, topographical features of the deposit were formed into a 3D geological model (Fig. 2). Finally, geostatistical studies were conducted and subsequently, the IDWAM combined with variography was proposed in order to build the RQD block model based on the dataset (raw data). The final block model was produced by applying an upper and lower filter using RockWorksTM software, based on the topographical surface data (Fig. 2) and borehole data collar heights. Those voxels located above the upper (based on topography) and below the lower and bed rock (based on collar heights) are not considered into the RQD block model construction [3,[42][43][44]. The Kahang deposit was modelled with 489,927 voxels and each voxel has a dimension of 4 m × 4 m × 10 m, corresponding the project dimensions of 600, 660 and 780 m, in the X, Y and Z directions, respectively. This has been done using statistical characteristics of the deposit geometrical particulars consisting of mean, median and median absolute deviation (MAD) [27,31,45,46]. For validation purposes between raw and estimated data, jackknife analysis was utilised.

Inverse Distance Weighted Anisotropic Method (IDWAM)
Inverse Distance (ID) is a common gridding and estimation method. Following this, the value assigned to a voxel is a weighted average of either all of the data points or a number of directionally distributed neighbours. The value of each of the data points is weighted with respect to the inverse of its distance from the voxel [17,30,31].
Inverse Distance Weighted Anisotropic method (IDWAM) is an interpolation of scattered points (RQD samples in this scenario) that estimates voxel values by averaging the values of sample data points in the neighbourhood of each processing voxel. IDWAM has a critical assumption that the interpolating surface is mostly influenced by the nearby points and less by the more distant points. The interpolating surface is a weighted average of the scatter points and the weight assigned to each scatter point diminishes as the distance from the interpolated point to the scattered point increases. The main advantage of the IDWAM is to produce a smooth and continuous grid and does not exaggerate extrapolations beyond the given data points [32]. Therefore, the IDWAM is recommended for geochemical mapping where the data boundaries (RQD populations) are critical for threshold-based target separation [19]. The range of RQD values will be smaller than the raw data range meaning that highest RQD values will be less than the maximum of raw data, and the lowest grade values will be greater than the minimum data point [47]. A general form of finding an interpolated value u at a given point x based on samples = ( ) for = 1,2, … , using IDW is an interpolating function: where x denotes an interpolated (arbitrary) point, xi is an interpolating (known) point, d is a given distance (metric operator) from the known point xi to the unknown point x, N is the total number of known points used in interpolation and p is a positive real number, called the power parameter (e.g., an exponent of "2" = Inverse Distance Squared, "3" = Inverse Distance Cubed). The greater the value of the exponent, the less influence distant control points will have on the assignment of the voxel value.
The disadvantages of conventional IDW methods are choice of weighting function which may introduce ambiguity especially where a fixed search radius requires a neighbourhood distance and a minimum or maximum number of points.

Statistical characteristics
Figure for the original data sets used for RQD values has been generated using MATLAB software, as depicted in Fig. 3. In this deposit, 14979 RQD samples have been measured from 42 boreholes carried out in the deposit (Fig. 3). RQD histogram provides a means for quickly evaluating the range of density and RQD values for a selected data set without creating a 3D solid model in order to illustrate the highest, lowest, sum, or average data values [6,9,31]. These histograms are used to read a single column of data (RQD or density) from a data set to determine the frequency or percentage of the total number of measurements for that variable/attribute that falls within each userdefined grouping. The RQD histogram is not normally distributed which follows a bimodal distribution, with average of 48% (Fig. 4).   (Table 1). In addition, with respect to the RQD histogram, the greatest frequency of the RQD data corresponds to the excellent rocks within the deposit which will result in the highest stability of the final pit slopе.

Block modelling
Selecting an accurate voxel size for a reserve/resource block model is important for minimising errors [30,47]. This issue has been investigated for estimated block model utilising various geostatistical methods such as ordinary kriging (OK) and IDW. Results achieved by the estimation methods are associated with the identification of voxel size in block modelling [27,48]. Utilising a larger voxel size will maximise the averaging effect in the estimated block model in terms of concentrations. Additionally, a smaller voxel size will show more details, but potentially more error in an anisotropic environment [31,49]. Furthermore, huge voxel size in the block model changes the higher or lower RQD values of RQD block model by smoothing of these points with high or low values within a large voxel. As a result, it is important to allocate an optimal voxel size considering the deposit geometry and drilling pattern (Figs 2 and 3) [50][51][52][53]. David (1970) proposed an applicable method for an operation based on geometrical particulars for different ore deposits and grid drilling [27]. Accordingly, voxel dimensions are determined as follows: Length and width of each voxel is equal to between half and quarter of the distance between the drill cores according to along the least variability deposit.
Height of each voxel is delineated due to the type of the deposit. In 'massive' deposits such as magmatic deposits (e.g., porphyry deposits), the parameter is equal to the height of excavating benches in the open pit mines [3,31].
The 2D map of 42 boreholes drilled in the Kahang deposit was constructed by RockWorks™ v.15 software (Fig. 5). As can be seen, the grid drilling pattern within this deposit is not homogeneous and systematic [17]. As a result, the Kahang deposit was modelled with 489,927 voxels; each voxel having a dimension of 4×4×10m in the X, Y and Z directions based on the geometrical properties of the deposit and grid drilling dimensions [27]. The RQD 3D model is evaluated by IDWAM, which can improve the interpolation of voxel values (RQDs) that lie between data point clusters and can be useful for modelling drill-hole based data in different types of deposit [51].

Figure 5. Location of boreholes sampled for density and RQD
Directional and non-directional searching in this method can improve the interpolation of voxel values that lie between data point clusters and be useful for modelling drill-hole based data in the stratiform and massive ore deposits. In this paper, a combination of IDWAM and variography has been used in order to generate a block model in terms of RQD values based on the following criteria: 1-The grid drilling pattern is irregular and nonsystematic (Figs 2, 3 and 5), with an especially high drilling density in the NE part of the deposit, and low density in the NW part (e.g., two isolated boreholes, as depicted in Fig. 5). Moreover, the grid drilling pattern has an anisotropic geometrical shape.
2-There are too many scattered drill holes in the marginal parts of the deposit which leads to a lack of data.
3-Trends of RQD values in X, Y and Z show that there is no association between ore grade and X-Y location or depth within the deposit (Fig. 6), indicating again that "Universal Kriging" is not appropriate for this deposit.

Application of IDWAM
For RQD, an experimental variogram in horizontal (with lag spacing of 15 m) and vertical (with lag spacing of 8 m) orientations was produced via MATLAB software with respect to log transformation of RQD values (raw data), as shown in Fig. 7. From this, the horizontal (Azimuth: 0 and Dip: 0) and vertical (Azimuth: 0 and Dip: -90) ranges for RQD are 120 m and 270 m, respectively. The theoretical variogram for RQD is as follow: (ℎ) = 135.71 + 190 ℎ (20,20,20) (3) +499.88 ℎ (120, 120, 270) 3D block model for RQD was generated by IDWAW using the RockWorksTM software package (Fig. 8). The horizontal and vertical ranges of RQD, 120 m and 270 m, respectively with respect to Equation 3 based on the conducted variograms (Fig. 7), were imported into the horizontal and vertical distance cut-offs in the "Solid Modelling Options" of the software. In order to achieve this, the following tasks in their relative order were carried out: The Weighting Exponent value was determined as being equal to 2 in order to prompt to enter a real number value for the Inverse-Distance exponent. Number of neighbouring points were defined between 3 and 15 data points that were to be used when computing the voxel value. The horizontal and vertical ranges (known as 'Cut-offs Distances' in the RockWorksTM software) were recognised based on the combined variograms with lags' spacing of 15 m and 8 m for horizontal and vertical directions ( Fig. 7 and Equation 3 with respect to the theoretical variograms). Following this, 3D block model for RQD was generated by IDWAM, as depicted in Fig. 8.
According to RQD block model in Kahang porphyry deposit, the RQD value for excellent rocks (RQD > 90) is situated in the central and NW parts of the deposit. Voxels including good (75 < RQD < 90) and fair (50 < RQD < 75) rocks are located in the central, eastern and NW parts of the deposit according to the RQD block model and classification of Rock Quality Designation, Deere and Miller rock classification (1966) [19] (Fig. 8 and Table 1).

RQD Spatial Distribution Models with Reliance on Deere and Miller Classification
The RQD block model (Fig. 8) is delineated to 5 populations with respect to Deere and Miller rock classification (1966; Table 1) for better understanding of RQD spatial distribution within the generated block model. To do this, the RQD block model was distinguished using a mathematical filter facility within the RockWorksTM software which is called "Boolean data type". As a result, the studied 3D model (Fig. 8) was allocated with binary codes (zero or one). Consequently, allocated zones (voxels) with RQD values with the code number of zero are removed and the zones with the code number of one will remain and appeared in the 3D models (Fig. 9). Based on the Fig. 9a, the majority of the deposit consists of very poor zones which is scattered within the Kahang deposit. As can be seen from Fig. 9b, poor zones are distributed within the most parts of Kahang deposit especially in the NE part. The fair and good zones with RQD values of 50-75 and 75-90, respectively, are present along NE-NW trend, as depicted in Figs 9c and 9d. Excellent zones in terms of RQD occur in the central and NW parts of the deposit (Fig. 9e).

Validation Processes
There are so many interdependent subjective decisions in a geostatistical study that it is a good practice to validate the results obtained by the estimation method (IDWAM in this scenario) prior to any production run. The generated block model is validated by re-estimating known values under implementation conditions, including the variogram model, estimation method and search strategy, as close as possible to those of the forthcoming production run [53][54][55]. For validation propose of the estimation, jackknife was used. The term jackknife applies to resampling without replacement, i.e., when alternative sets of data values are re-estimated from other non-overlapping data sets. The jackknife analysis in the Kahang deposit indicates that the correlation coefficient between original data and RQD estimated is +0.78, as shown in Fig. 10. In this figure, the diagonal of the square plot (black line) and the linear regression (blue line) was derived and calculated using RockWorksTM program. The jackknife analysis in the Kahang deposit for RQD data indicates that the correlation between original data and estimated RQD via IDWAM is over +0.7. In this figure, the diagonal of the square plot (black line) and the linear regression (blue line) has a small angle and difference. However, there are some numbers of unexplained variances.

Conclusions
Results achieved by the proposed IDWAM estimation with combined variogram indicated a good correlation equal to 78% between RQD estimated and raw data. Therefore, it may be resulted that the proposed method can be applied for interpolation of geomechanical data in different types of ore deposits while dealing with nonuniform grid drilling like what it is in Kahang deposit. The RQD block model obtained by the proposed estimator illustrate that the excellent zones considering RQDs higher than 90% occur in the central and NW parts of the deposit. In addition, the good zones with RQD values between 75-90% are present in the central, NE and NW parts of the studied area. It is believed that a final pit slope geometry and ultimate pit limit depend on the geomechanical properties such as RQD. However, an awareness of the spatial variability of parameters such as RQD can be used to investigate geotechnical characteristics. This can then be utilised to assess potential slope stability and be incorporated into a geotechnical risk model for the final pit geometry. Regions of high RQD may be targeted as offering greater potential for increased slope angles or locations for siting of critical haul roads. Regions of lower RQD should, be avoided for final pit limits. From a slope stability point of view, it may be expected that anyone examining RQD block model would consider at least multiple domains for slope stability assessments.
However, this is highly recommended that more geotechnical characterisations are necessary to be examined for any potential influence of the 3D fracture network and any major discontinuity-controlled instability.