Stream Sediment Geochemical Mapping Using Kriging Methodwith Optimized Parameters(Case

Stream sediment geochemical mapping using kriging methodwith optimized parameters(case study: Khoy district, northwest of Iran)

Farhad M. Torab1 and Sara kasmaee2

1- Assistant Professor, Department of Mining and Metallurgical Engineering, Yazd University, Yazd, Iran

Corresponding author: e-mail: , tel: +98-351-8122521

2- Department of Mining and Metallurgical Engineering, Yazd University, Yazd, Iran e-mail: ,tel: +98-351-8122521

Abstract:

The aim of geochemical mapping in this paper is to obtain regional background and delineate geochemical anomalies for mineral exploration using the combined application of data analysis and optimized kriging. Multi-element geochemical dataset from stream sediments in Khoy district (northwest of Iran) was used to illustrate this application. The Factor analysis was first utilized to establish the different associations among the variables resulting 4 distinct factors and afterwards to obtain a homogeneous group for applying the geostatistical methodology. Taking advantage of high correlation within the group of elements and the components, we used PC1 and PC3 as new variables for numerically identifying potential areas. The first group composed of Cu, Fe, U, Hf, V, Sn,Ta and Th,strongly correlated together,were selected as first principal component (PC1).The third group consists of Cr, Ni, Co and Mg represents third principal component (PC3). The geostatistical approach involved computing variograms of PC1, PC3 and subsequently mapping of this variables using ordinary kriging.Optimization of the kriging parameters while mapping the geochemicalanomaliesis a suitable strategy for defining themost precise estimation.The result of kriging efficiency analysis demonstrates its value as a novel method to utilize the best estimating kriging parameters in geochemistry.When attempting to differentiate between anomalous samples and those belonging to geochemical background, choosing the best cell size and appropriate kriging parameters are essential. By optimization of kriging parameters, the best cell size 1000 m is chosen in order to increase accuracy, indicating better delineation of geochemical anomalies.To sum up, due to the optimization results,the 74.67% kriging efficiency was attained for PC1 and 67.34% kriging efficiency for PC3 in geochemical mapping.

Key words: geochemical mapping, ordinary kriging, kriging optimization, kriging efficiency, cell size

Introduction:

Regional geochemical mapping is the study of the chemistry of the surface of the earth. It can be considered as an important tool in environmental studies, geological mapping and mineral exploration. Objectives of regional geochemical mapping include 1) elucidation of fundamental earth system processes, 2) demonstrating geochemical anomalies for exploration of new mineral resources, 3) determination of geochemical background values, and 4) identification of natural or manmade chemical contamination (Zumlot et al 2009). Geochemical mapping requires a good understanding of geochemical variables and proper use of mathematical and mapping techniques (Chaosheng Zhang et al 2013). The multivariate and regionalized character of geochemical variables makes them an interesting candidate for numerical analysis using both geostatistics (Matheron, 1970; Journel and Huijbregts, 1978) and data analysis methods (Benzecri, 1973; Davis, 1987) in order to identify geochemical mapping. In practice, it is difficult to determine the suitable cell size in geochemical mapping. The effects of cell size on the structure of the geochemical mapping and the determination of the thresholds are seldom addressed (Zuo 2012). The selection of cell size is usually based on the density of sampling or mapping scale, or the structure of point pattern. The interpolation with different cell sizes may result in different concentration frequency distributions and different textures of geochemical patterns (Hengl, 2006).

For geochemical mapping multivariate statistics can be reasonably used when various elements should be considered. In this sense, the results become more important than isolated anomalies of individual elements (Jimenez-Espinosa et al 1992). Recently different methodologies have applied as interpolation technique and anomaly discrimination such as: moving weighted median (MWM), kriging, inverse distance weighted (MIDW), and concentration–area (C-A) fractal method (Lima et al, 2008).

This paper aims to interpret and deals with the concentrations of 52 elements obtained through chemical analyses of stream sediments in Khoy district located in northwest of Iran. Representative geochemical maps of the Khoy region are presented, and the lateral variations of elements are discussed with respect to the geochemical characteristics of the ophiolite complex rocks, distribution of mineral deposits and igneous activities.

In this research, by applying ordinary kriging analysis, it is possible to identify geochemical anomalies of different components, which are representative of a related group of elements. Thedifferent parameters used in an estimate can be studied and optimized to improve accuracy of geochemical mapping. For the first time, the effects of kriging parameters and cell size are explored in geochemical mapping and by optimizationof the parameters,the optimum kriging efficiency is calculated in estimation of each principal components.

2. Case study presenting: geology and sampling

2.1. Geology

Studies of regional geochemical mapping using stream sediments have been conducted in northwest of Iran; Khoy district. The Khoy Mountains are spread with an NE–WS trend, the highest peak is in southwest, Ervin mountain (3622 m), and the lowest part is the Khoy basin near the Khoy city (1080m), (Fig. 1). The Main rivers; Aghchai and Aland flow eastward and northward along the Ervin Mountains, alternatively. On the south side of the Ervin Mountains, rivers flow northward to the Aras river. The geological survey of Khoy has undertaken a geochemical mapping, aiming at documenting the elemental distribution in surface areas and providing a database for geochemical assessment. Stream sediments with a sampling density of one sample per2.7 km2 were taken with sample analyses conducted for more than 50 elements. These studies have indicated that chemical variations in the stream sediments are strongly controlled by geology and the distribution of mineral deposits. The underlying lithology of the study region comprises of quaternary sediments, and in some parts of the area there are conglomerate with intercalation of sandstone, Pliocene intrusive rocks (quartz monzodiorite, granodiorite with minor diorite, pegmatitic dikes), rocks of ophiolitecomplex (Mesozoic basaltic lava flows, basaltic pillow lava and partly serpentinized ultramafic rocks),amphibolite and metamorphic rocks with unknown age.According to the complexity of geology in this area, it is essential to use mapping of geochemical anomalies to prospect the potential areas for further studies.

2.2. Sampling, sample preparation and analysis

Throughout the exploration project using stream sediments, 844 samples were collected and analyzed for geochemical mapping (Fig. 2). The grain size is varied, but sands and gravels are typically abundant at the sampling sites. Samples were dried in air and sieved through a 0.18 mm (80 mesh) stainless steel screen. After this, they were crushed, homogenized and sieved, retaining the <200 mesh fraction for chemical analysis. Analyses were performed by atomic absorption spectroscopy (AAS) for Au, and an inductive coupled plasma atomic emission spectrometer (ICP-AES) for other elements.

To check the quality of the analysis, duplicate samples were analyzed and the analytical errors as relative percentage difference (% RPD) were calculated as follows which arepresented in Figures 3.

% RPD = (S-D)/((S+D)/2)×100 (1)

Where: S = determinate value of the samples, and D = value of the duplicates (Guillén et al 2011). Nearly 5% of the samples were analyzed in duplicates using internal controls.

The accuracy and the precision errors (RPD: relative percentage difference) were estimated as a filter choosing the appropriate elements. The errors less than 30% denote as an acceptable confidence in the results.

Statistical analysis

The objective of this work is to describe the geochemical characteristics and portray the spatial distribution and variation of elements in stream sediments.Grunsky (2010) described a systematic approach to evaluating geochemical data that involves the examination of geochemical data as both individual elements and multivariate associations. For interpreting data, all acceptable stream sediment data –have chosen in part 2.2-are used to interpret the high statistics, box-plots, histograms, pearson correlation analysis, cluster analysis and robust principal factor analysis (Reimann et al., 2008; Grunsky and Kjarsgaard, 2008; Grunsky, 2010). The purpose of using these techniques was to investigate the structure, distributions, trends and element associations within the data set, thus providing insights into the underlying geological, physical and geochemical processes that are important in controlling the stream sediment geochemistry (Lapworth et al,2012). In the first step, univariate statistical analysis was used to investigate the outliers and then modifying them using Humpletest (Linsingeret al.,1998; Doerffel,1962). Secondly, plotting histograms showed data distributions due to necessity in symmetric shape in data distributions and comparable magnitude value range for multivariate analysis including cluster analysis (Reimann et al., 2008).Figures 4 and 5 show two examples,the histograms for raw data samples and then the histogram of normalized data after modifying the outliers.

Significant differences in range value could cause spurious pattern due to dominance of particular components. Thus, when data are not following the condition above it is important to apply some process such as data transformation before applying multivariate analysis (Nugraha et al,2011).The geochemistry data was stored in a relational MS EXCEL sheet and statistical analysis and plotting was carried out using SPSS software package. The summary statistics for the stream-sediment geochemistry after modifying the outliers and normalizing someelements by Log-normal transformation are shown in Table 1.

Clustering techniques were used to separate the individuals into natural groups; by means of proximity criteria, the discrepancies between the individuals in a single group were minimized while the distances separating the groups from each other were maximized (Martinez et al., 2007). In order to discriminate different groups of elements as tracers of a natural source, an explorative hierarchical cluster analysis was performed on the available data set. The results obtained (Fig. 6) enabled four main groups based on a visual observation of the dendrogram, namely C1-C4according to their respective predominant lithology.

Principal Component Analysis (PCA) has been widely used to explore datasets that are large and spatially diverse which otherwise makes the determination of spatial characteristics difficult (Reid and Spencer, 2009).When this technique is applied to a geochemical data set, it is possible to obtain several components, as linear functions of the original chemical elements, which identify subsets of correlated variables. Some of these components can be used for studying a specific group of variables, giving conclusions about an association of elements, which is geochemically more significant than the study of individual variables. In this way, we can use one of the principal components for identifying geochemical anomalies by means of a geostatistical approach (Jimenez-Espinosa et al 1993).To complement geochemical characterization, there have been taken into account the data resulting from statistical analysis together with the spatial distribution of elements obtained through the compilation of baseline geochemical maps. In this study a Varimax rotation was included to maximize the loading of variance onto each principle component to aid data interpretation (Kaiser, 1958; Zhou et al., 2004).Table 3 gives the results of the PCA factor loading for the Khoy stream sediments data. The extracted components PC-1, PC-2, PC-3 and PC-4 represent 32%, 25%, 8% and 5% of the total variance of the geochemical data, respectively.

The 1st component (PC-1) explains the dispersion of correlated Cu±Fe±U±Hf±V±Sn±Ta±Th and their association with major igneous rocks of the area. The Be±Sm±Ce±K±La±Nd±P combination appears in the PC-2, highlights mostly the REE elements in the study area. The third component (PC-3) extracted from the stream sediment geochemical analysis represents the Co±Cr±Mg±Ni associations indicating the ultramafic rocks and ophiolite complex. The last component (PC-4) consists of Al±Dy±Nadescibes the lithological backgrounds with less importance in the area. According to the prospecting studies carried out in Khoy district and geology maps of the area, components of PC1 and PC3 have been chosen as the most effective factors in determining the geochemistry potential areas.

The 1st component (PC-1) explains the dispersion of correlated Cu±Fe±U±Hf±V±Sn±Ta±Th and their association with major igneous rocks of the area. The Be±Sm±Ce±K±La±Nd±P combination appears in the PC-2, highlightsmostly the Rare Earth Elements(REE) in the study area. The third component (PC-3) extracted from the stream sediment geochemical analysis represents the Co±Cr±Mg±Ni associations indicating the ultramafic rocks and ophiolite complex. The last component (PC-4) consists of Al±Dy±Nadescibes the lithological backgrounds with less importance in the area. According to the prospecting studies carried out in Khoy district and geology maps of the area, components of PC1 and PC3 have been chosen as the most effective factors in determining the geochemistry potential areas.

4. Geostatistics

Geostatistics provides an advanced methodology which facilitates quantification of the spatial features of geochemical sedimentary samples and enables spatial interpolation. Geostatistics is based on the theory of a regionalized variable (Matheron, 1963), which is distributed in space (with spatial coordinates) and shows spatial autocorrelation such that samples close together in space are more alike than those that are further apart. Geostatistics uses the technique of variogram (or semi-variogram) to measure the spatial variability of a regionalized variable, and provides the input parameters for the spatial interpolation of kriging (Krige, 1951; Webster and Oliver, 2001).

It relates the semi-variance, half the expected squared difference between paired data values Z(x) and Z(x+h), to the lag distance h, by which locations are separated:

(2)

Where Z(xi) is the value of the variable Z at location of xi, and N(h) is the number of pairs of sample points separated by the lag distance h. For irregular sampling, it is rare for the distance between the sample pairs to be exactly equal to h, therefore, the lag distance h is often represented by a distance band.

A variogram plot can be acquired by calculating variogram at different lags. Then, the variogram plot is fitted with a theoretical model, such as spherical, exponential, or Gaussian models. The fitted model provides information about the spatial structure as well as the input parameters for kriging interpolation (McGratha et al, 2004).

Experimental variograms of PC1 and PC3 were calculated in several different spatial directions in order to analyze presence of anisotropies. Table4 presents the fitted spherical variogram model parameters to be used in ordinary kriging modeling in two different factors. As it is demonstrated in table 4, the two PCs show different variogram parameters especially anisotropy ratio, which these differences will be erroneously ignored in conventional estimations in routine geochemical mapping.

Kriging (Krige, 1951) is regarded as the best linear unbiased estimator, which is a process of a theoretical weighted moving average:

(3)

where is the value to be estimated at the location of Xi, Z(xi) is the known value at the sampling site xi. Altogether there are n samples used for the estimation, and the number n is selected based on the size and shape of the moving ellipsoid window and depending on the range and anisotropy ratio selected as user definition. Contrary to other methods (such as inverse distance weighted), the weighting function γi is no longer arbitrary, and it is calculated based on the parameters of the variogram model. To ensure that the estimate is unbiased, the weights need to sum to one:

(4)

The estimation errors (or kriging variances) need to be minimized. As in conventional statistics, a normal distribution for the variable under study is desirable in linear geostatistics (Clark and Harper, 2000). Even though normality may not be strictly required, serious violation of normality, such as too high skewness and outliers, can impair the variogram structure and the kriging results. However, in this case, the parameters PC1 and PC3, as the important factors, are essentially normal standard.

Krige (1996) presents a practical analysis of the effects of spatial continuity and the available data within the search ellipse as it affects measures of conditional bias. The two parameters he suggests using to investigate whether the block size used for grade estimation is appropriate are kriging efficiency (KE%) and regression slope (R) which can also be used to calibrate the confidence in block estimates and are given as follows:

KE = (BV-KV)/BV (5)

R = BV-KV+ |μ|/ BV – KV + |2μ| (6)

Where:

BV = theoretical variance of blocks within domain;

KV = variance between kriged grade and true (unknown) grade, ie kriging variance;

μ = Lagrange multiplier

Perfect estimation would give values of KV = 0, KE = 100% and R=1 (Snowden 2001).

Many methods of classification based on kriging variances, kriging efficiency or regression slope have been observed. Not all of them are valid, as pointed out by (Glacken 1996) and it may be worth noting a few poor practices (Snowden 2001). In this study,we have considered the kriging efficiency as the best function to choose the most suitable estimation parameters in geochemical mapping.

Selection of optimized parameter in geochemical mapping

One of the most important stages in mineral exploration and environmental studies is geochemical mapping. Whereas geochemical point data are collected from stream sediments, finding undiscovered mineral deposits and monitoring environmental changes is a fundamental step (Zuo 2012). Hence, the cell size should be fine enough not only to indicate the potential areas but also shows an appropriate interpolation result. In practice, the appropriate cell size can be determined from the density of samples and mapping scale, or the structure of pattern of points (Zuo 2012), estimation method and its input parameters. In the first step, the study area was divided into various grid systems with the cell size of 100 m, 400m, 600m, 800m, 1000m, 1200 and 1400m. Two distinct maps showing the distribution of geochemical halos and anomalies were obtained for PC1 and PC3 as the main factors by applying the ordinary kriging method. The Ordinary kriging interpolation and mapping were applied for 7 block sizes, tested with various numbers of samples and different search distances. According to the previous studies on investigating the appropriate block size using the kriging efficiency (KE%)- above 50 percent is advisable by Krige- ; it has chosen here to calibrate the confidence in optimization of kriging parameters and even choosing the precise geochemical maps. Therefore, all the tested estimations are compared by the results of their mean kriging efficiency (Fig 10and 13). The regularization of OK results shows that not only the cell size influences on kriging efficiency but also it might be so sensitive to kriging parameters.