Sequential Gaussian Simulation in the Sungun Cu Porphyry Deposit and Comparing the Stationary Reproduction with Ordinary Kriging

Main aim of this study is stationary (first and second order) reproduction assessment of sequential Gaussian simulation (SGS) compared to ordinary kriging (OK) in the hypogene zone of the Sungun Cu porphyry deposit (NW Iran) based on drillcore data utilizing both SGS and OK. In comparison with OK, geostatistical stochastic simulations are able to assess spatial and uncertain distribution of desired variables. In simulation by SGS method, data were transferred to standard normal and then simulated 10 times (10 realizations were created in this scenario. All of the realizations were honor to histogram and variogram of samples, so all realizations are valid but derived cumulative distribution functions from E-type map and kriging could not reproduce the stationary well.


Introduction
For better understanding of the grade variability, it is very important to predict the spatial distribution of the grade within the desired ore deposit [1]. Geostatistics has been used for spatial variability characterization and prediction of grade over the three last decades. Ordinary kriging (OK) is the most useful geostatistical estimation technique which is called also the "best linear unbiased estimator" [2,3,4]. The most important negative characteristics of moving averages estimators such as kriging is smoothing effect and reducing the range of variation of the variables. Geostatistical simulation is widely used to overcome this problem and avoiding the smoothing effect of such estimation methods [1].
Geostatistical stochastic simulation could also assess local and spatial uncertainty of the estimates [2,3]. Conditional stochastic simulation is designed initially to overcome the smoothing effect of kriging estimator especially when mapping sharp or extreme spatial discontinuities is looked for [3,4]. The simulation algorithms take into account both the spatial variation of actual data at sampled locations and the variation of estimates at unsampled locations [5]. It means that stochastic simulation reproduces the sample statistics (histogram and semi-variogram model) and honors sample data at their original locations. Therefore a stochastic simulation map represents the spatial distribution of the attribute more realistic than a kriged map [6,7].
Many stochastic simulation algorithms have been proposed and sequential gaussian simulation (SGS) among them is widely used due to being fast and straightforward in reconstructing Conditional Cumulative Distribution Function (ccdf) [8,9,10,11].
In this paper, Sequential Gaussian Simulation (SGS) and Ordinary Kriging were utilized for assessing the stationary reproduction ability of each simulation/estimation methods based on Cu values in the hypogene zone of Sungun porphyry copper deposit (PCD), NW Iran.

General Specifications of Study Area
In Iran, all known porphyry copper mineralization occurs in the Cenozoic Sahand-Bazmanorogenic belt (Fig. 1A). This belt was formed by subduction of the Arabian plate beneath central Iran during the Alpine orogeny and hosts two major porphyry Cu deposits. Sar-Cheshmeh, Sungun, and a number of sub economic porphyry copper deposits are all associated with mid to late Miocene diorite-granodiorite to quartz monzonite stocks [12,13,14]. The deposit is situated in the Cenozoic Urumieh-Dokhtar magmatic belt of the Zagros orogen extending from NW to SE Iran [15].

Geologic Setting of the Sungun Deposit
The Sungun PCD is hosted by a diorite-granodiorite to quartz monzonite stock, located 75 km northwest of Ahar in the Azarbaijan province of northwestern Iran (Fig. 1A) [13]. The stock is part of the Sahand-Bazman igneous and metallogenic belt (northern Iran), a deeply eroded Tertiary volcanic field, roughly 100 by 1,700 km in extent (from Turkey to Baluchistan in southern Iran), consisting mainly of rhyolite and andesite, with numerous felsic intrusions. The volcanics were laid down unconformably over folded and eroded Upper Cretaceous andesitic volcanic and sedimentary rocks (~500 m thick)

Sungun Stock and Peripheral Intrusive Rocks
The Sungun stock is a complex intrusive body, which crops out over an area of about 1.5 by 2.3 km. The stock consists of three different intrusive phases: (1) monzonite-quartz monzonite, (2) diorite-granodiorite, and (3) andesite and related dikes, listed in order of emplacement [13,16]. The monzonite-quartz monzonite is volumetrically the most important and includes most of the western part of the intrusive complex at the current erosional surface. The diorite-granodiorite is volumetrically the next most important and hosts most of the mineralization. These two intrusive phases are cut by monzonitic and andesitic dikes, which in the northern and eastern parts of the Sungun stock are also locally mineralized [17] The dykes trending NNW-SSE, dip steeply to the west and have thickness from a few centimeters to several tens of meters (Fig. 1B). What is important and problematic is that the dykes appear to have acted as a barrier to hydrothermal and supergenes processes [18]. Two distinct supergene alteration zones are recognized at Sungun (1) oxidized and leached zone and (2) supergene sulfides zone [17].

Alteration and Mineralization
Hydrothermal alteration and mineralization at Sungun PCD are centered on the stock and were broadly synchronous with its emplacement. Early hydrothermal alteration was dominantly potassic and propylitic and was followed by later phyllic, silicic, and argillic alteration [16,18].

Hypogene Zone
Hypogene copper mineralization was introduced during potassic alteration and to a more extent during phyllic alteration [18]. During potassic alteration, the copper mineralization was deposited as chalcopyrite and minor bornite; later hypogene copper mineralization deposited mainly chalcopyrite. Alteration of feldspars and biotite (from potassically altered rocks) was accompanied by an increase in sulfide content outward from the central part of the stock. The maximum Cu grade is associated with biotite, orthoclase, and sericite (potassic zone) while the pyrite content is highest (3-10 vol % of the rock) in the marginal quartz-sericite (phyllic) zone [16].

Sequential Gaussian Simulation (SGS)
Sequential simulation is a stochastic modeling algorithm that obtains multiple realizations based on the same input data [10,19]. This data could be either continuous or categorical. Regarding the data type, sequential indicator simulation, sequential Gaussian simulation (SGS) or direct sequential simulation will be used [20].
The most straightforward algorithm for generating realizations of a multivariate Gaussian field is provided by the sequential principle [4,8]. SGS demands standard Gaussian data with zero mean and unit variance, so for SGS, data are transformed to be Gaussian through a quantile transform [1]. Each variable is simulated sequentially according to its normal ccdf through a Simple Kriging estimation system. The conditioning data consist of all original data and all previously simulated values found within a neighborhood of the location being simulated [4,8].
The conditional simulation of a continuous variable z(u) in a Gaussian space proceeds as follows [ [3].
Regarding to a transformation to gaussian and then backtransform to original unit, statistical fluctuations are inherent in simulation but the fluctuations should be reasonable and unbiased in the mean and variance [4,8].
Following checks should be performed after having all nodes simulated: reproduction of (1) data values at data locations, (2) the original histogram, (3) the original summary statistics, and (4) the input covariance model [4].

Descriptive Statistics and Spatial Autocorrelation Analysis
Descriptive statistics and the histogram of copper grades from 19605 samples taken from 135 drill holes in the hypogene zone of the Sungun PCD (Fig. 2) are presented in Table 1 and Fig. 3, respectively Which show the distribution of Cu data is not normal with the Cu mean value of 0.445% (Fig 3A).   The Cu data have been transformed by using a normal score transformation and the statistics of transformed data (i.e. the mean value of 0 and the variance of about 1) check the correctness of the transformation.
The display of the histogram of new Gaussian variable also checks that the distribution is symmetric with minimum and maximum values of −3.88 and 3.88 respectively (Fig. 3B). Geostatistical studies and visualizations were done with SGeMS and Datamine studio. The experimental semi-variogram and a spherical model fitted to the raw and normal transformed data are presented in Fig. 4. The isotropic semi-variogram of raw data follows a spherical structure with a nugget effect of 0.038 (%) 2 which reaches to a sill of 0.22 (%) 2 at a range of 64.7m, as shown in Fig. 4A. The experimental semi-variogram of the normal scores (Fig. 4B) follows again a spherical model with a nugget effect of 0.138(%) 2 with range of spatial correlation is about 71m, which is similar to the one obtained for the raw data. That means the intrinsic spatial character of the Cu data does not vary with the normal transformation of data. The sill of the semi-variogram for the transformed data reaches to unity, as it should be to fulfill the second-order stationary assumption.
Anisotropy has also been investigated and modeled based on calculating the experimental semi-variograms of Cu and normalized value among horizontal directions with 30° angular increments and ±15° angular tolerance and one vertical direction. The results show a mild anisotropy in the azimuth of 150 (major axis) for normal values, as depicted in Fig.5A & B.

Conditional simulation
Based on SGS modeling, ten realizations of Cu spatial distributions are generated on a 25×25×25 (m 3 ) grid within the hypogene zone. Simulation is performed using the simple kriging estimator, and the semi-variogram model of Cu normal scores. Horizontal plan of four randomly selected realizations consisting of simulations number (sim 2, 5, 7 and 10) are displayed in Fig. 6. Each realization represents a realistic spatial distribution of Cu without a smoothing effect.
Four randomly selected realizations are checked to examine the sample statistics reproduction. The CDF's of all realizations of Cu distribution and also E-type and kriged map are displayed, as depicted in Fig. 7. Comparing these frequencies with the sample Cumulative frequency (solid red line in the same shape) reveals that the realizations reproduce the sample histogram, reasonably well (Fig. 7). The reproduction of the Cu raw data semi-variogram model by selected realizations is also proper.
Some acceptable discrepancies between different realizations and sample model called ergodic fluctuations which may have different reasons such as the algorithm used for the simulation, the semi-variogram model parameters and the amount of conditioning data to be utilized for the simulation [2].  In the case of SGS algorithm, the histogram and semi-variogram model reproduced over a number of realizations should be, on average, equal to the sample statistics (Fig. 8) [21].
Each realization well preserves the range of variation of the measured Cu data compared to OK map which is illustrated in Table 2. This reveals the smoothing effect, a typical property of kriging. This is also evident from the kriging variance, which is much smaller than the actual variance. OK does not reflect the true variability especially for high values (Fig. 9), and hence is not appropriate for the underlying outlook to this study. Summary statistics of E-type estimate map are however similar to those for the kriged map. On the other hand, results obtained by E-type have good correlation with estimation values derived via ordinary kriging.

Conclusion
The most important properties of Geostatistical Simulations are, producing a group of images which shows a range of possible events, calculating probable percentage of happening and also determining the risk in each step of process. The three dimensional deposit models were generated for Sungun PCD using SGS algorithm. All realizations could reflect the first and second stationary of Cu i.e. 10 realizations could reproduce the histogram and variogram of raw data that indicates the simulation process is validate. Up to here, all realizations are valid and all honour the first and the second order stationary of the data, but for mine-planning purpose, considering the calculated NPV of all realizations a ranking should be done through which the first rank would be the best realization for mine-planning purpose.