Geochemical and Behavioral Modeling of Phosphorus and Sulfur as Deleterious Elements of Iron Ore to Be Used in Geometallurgical Studies, Sheytoor Iron Ore, Iran

Sheytoor Iron Ore deposit is located in Yazd province of Iran (Bafq). The most abundant ore is magnetite, which can be seen in the form of mass and granular tissue in various forms of self-shaped, semi-self-shaped and amorphous. The main purpose of this study is to identify the geochemical relationship of phosphorus and sulfur elements and also three-dimensional modeling of mineralization of these elements in iron ore. In order to achieve the research goal, methods such as k-mean clustering technique, concentration-volume fractal as well as block modeling with kriging estimator and Inverse Distance Weighting (IDW) interpolator were used. The model of geochemical behavior of phosphorus and sulfur elements compared to iron is of great impor-tance because these two elements are known as deleterious elements in mineral processing and steelmaking processes, which are the post-mining stages. Existence of geochemical model and identification of elements’ behavior towards each other play a key role in optimizing mining operations in order to achieve geometallurgical goals. The results of this study are the three-dimen-sional model of mineralization of iron, phosphorus and sulfur elements, separation of phosphorus and sulfur mineralization communities and also pre-senting the model of enrichment community of these two elements. All the results are in line with geometallurgical studies and can optimize the next steps by optimizing the mining process.


Introduction
Iron is the most abundant metallic element in the earth's crust after aluminum.
The average abundance (Clark) of iron in the earth's crust is 4.56%. Ten times the enrichment of iron compared to Clark can form an economic deposit [1] [2].
The properties of iron have made it more widely known as the metal that gives us tools, and its ability to make it a major metal in technology using alloys and heat treatment suitable for any application. Iron is the most common metal in everyday life, always in the form of manufactured objects, and is usually covered with a protective coating. Concrete structures contain iron, and electrical machines, including transformers, are iron dependent [3] [4]. Cars are mostly made of iron, and many iron tools are covered with a protective layer of tin [5]. Fasteners, such as nails and screws used, are usually made of iron [6]. The deleterious elements of iron ore are phosphorus, sulfur and silica. Phosphorus causes steel to break. Must be removed or reduced [7]. The sulfur element is a harmful to steelmaking and softens steel and should be removed entirely. Silica is the most expensive annoying element of iron ore, which must be poured out of the melting pot in the form of melt and slag along with a large amount of iron [8].
Iron ore mineralization has existed throughout the geological period from the Proterozoic to the Quaternary, but the major sources of iron ore were formed during the Archaeozoic. In Iran, most sources of iron ore are formed in Bafgh-Posht-e-Badam zone, Sanandaj-Sirjan zone, volcanic-plutonic belt of North Khorasan (Sangan) and other tectonic-magmatic zones [9]. Geochemical behavior is considered as an important tool in identifying the relationships between elements and how deposits are formed [10]. Clustering methods can be of great help in the process of studying geochemical behavior. One of the most important clustering methods is the k-means technique. In this method, bulk data is divided into representative clusters, the members of each of these clusters have similar characteristics to each other [11] [12]. Using the K-mean clustering method, the concentrations of the elements in relation to each other can be predicted. Geochemical behavior of elements to each other can also be identified [13] [14]. Statistics is a powerful tool in geochemical studies that plays an important role in data modeling [15] [16]. Block modeling, in fact, provides a three-dimensional model using data on the concentration of elements [17]. In order to build a three-dimensional model of mineralization, based on the data obtained from drilling exploratory boreholes, various interpolation methods are used. Among the most widely used of these methods are Inverse Distance Weighting (IDW) and Kriging [18] [19]. Mineralization of an element involves diverse mineralization communities. These communities include communities background, anomalies, and enrichment. In two-dimensional (surface) studies, these communities can be identified by tools such as airborne geophysics and remote sensing (satellite image analysis). Each of these methods has already been used in several studies by researchers [20]- [25]. In-depth studies, we will also be able to separate these communities by modeling the concentration of elements  [34]. The results of this study can be the basis of geometallurgical studies in Sheytoor iron ore mine. Using these results, the mining process can be optimized. In this way, the production of the mine can be commensurate with the feed of the processing plant and the steelmaking process.

Geolocation of Sheytoor Iron Deposit
The Sheytoor iron deposit is in Yazd province of Iran and is located 78 km from Bafq city, near the gazestan village. The area of the deposit is approximately 28 km 2 . You can see the eastern view of the Sheytoor region in Figure 1 [35]. Due to the geological phenomena in this region, this region can be considered as one of the geotourism attractions [36].

Geological Setting
The area is located in Yazd province and in central tectonic zone of Iran, which is a small part of the Poshtebadam-Bafq block. In terms of lithology, it has a wide range of alkaline (gabbro-diabaz) to acidic (rhyolitic and rhyodactic) rocks, and includes intrusive, semi-deep and extrusive types. Eastern view of the region and the position of the drilling equipments is shown in Figure 2.
In Figure 2, the locating of 1:1000 map of Sheytoor deposit and 1:20,000 map of the gazestan in the topographic map of 1:50,000 Sheytoor is shown [35] [37].

Stratigraphic Units
In terms of lithology, the area consists of the late precambrian to cambrian volcanic-sedimentary rocks and igneous rocks after that. In the central parts of the

Mineralization
Mineralization of iron-phosphate is closely related to the region's rocks including gabbro, monzo gabbro, diorite, syenite and diabase which is seen in dark green color. The main ore deposit is as apatite-magnetite mineralization.

Raw Data Preparation Methods
Before using raw data must identified and replaced censored and outlier data.
Censored data is said to be data that among them, due to the high sensitivity limit of measuring devices, a number of data are found to be smaller than the device sensitivity limit. Such data can make statistical problems, because, firstly, statistical methods require a complete set of non-censored data, and secondly, in some cases, such as anomaly separation from background and relative  measurements, the existence of censored data leads to inappropriate evaluations.
If the censored data are identified and replaced, the amount of background and intensity of the anomalies will be calculated more accurately [38].
As the existence of censored data among geochemical data leads to errors, outlier data also have the same effect on the results. In statistics, an outlier is an observation point that is distant from other observations [39] [40]. An outlier may be due to variability in the measurement or it may indicate experimental error; the latter are sometimes excluded from the data set [41]. An outlier can cause serious problems in statistical analyses. Outliers can occur by chance in any distribution, but they often indicate either measurement error or that the population has a heavy-tailed distribution.
Several methods can detect and replace censored and outlier data. In this study, a simple method is used to replace sensor data. In this method, the values of less than sensitivity limit are replaced by 3/4 of data value. The main problem of this method is that it is by no means influenced by the statistical parameters of the data society and is merely a function of the sensitivity limit of the measurement method [38] [42].
In order to identify and replace the outlier data, the doerffel method was used. Doerffel has prepared a graph for determining the threshold of outlier data values, which is provided for two levels of significance of 5% and 1% (see Figure 4) [35] [42].
To perform the doerffel test, the average ( x ) and standard deviation of the data (s) is calculated regardless of the largest amount of data. Then the largest amount of data (x A ) is considered to be outside of the row if it is true in the following Equation (1): where "g" is the threshold limit for outlier values which can be calculated by graph shown in Figure 4.

Estimation Methods in Modeling
In this study, in order to grade estimation, methods kriging and inverse distance weighting (IDW) were used.   Kriging: The basic idea of kriging is to predict the value of a function at a given point by computing a weighted average of the known values of the function in the neighborhood of the point. In geostatistical models, sampled data is interpreted as the result of a random process. The fact that these models incorporate uncertainty in their conceptualization doesn't mean that the phenomenon-the forest, the aquifer, the mineral deposit-has resulted from a random process, but rather it allows one to build a methodological basis for the spatial inference of quantities in unobserved locations, and to quantify the uncertainty associated with the estimator. The kriging function is defined as follows (2): λ i is the weight associated with the value of the x-variable x at the point i in the case where Σλ = 1 [43].  Inverse Distance Weighting (IDW): Inverse distance weighting (IDW) is a type of deterministic method for multivariate interpolation with a known scattered set of points. The assigned values to unknown points are calculated with a weighted average of the values available at the known points. IDW is an interpolation method in which estimation is performed based on the values of the nearest points to the point of the weighted inverted distance. The calculation method of the IDW is as follows Equations (3) and (4) [44]:

Three-Dimensional Modeling
3D modeling is one of the most common types of modeling to understand the results in a form of ore below the surface of the earth. Considering the threshold limit, it is possible to see the separated part of the deposit, which is considered as an anomaly. The most common type of 3D model for the description of the ore and the implementation design is "Regular 3D Fixed Block Model". Figure 5 shows an example of a block model [45].

K-Means Clustering Method (Behavioral Study)
K-Means clustering is a method for classifying voluminous data. In this method, using clustering algorithm, data is divided into optimal classes. In this study, using the K-means method, the geochemical behaviors of iron interfering elements are also identified.
The k-means algorithm starts with a given value for K (number of classes) and tries to estimate the following cases:  Finding the points as centers of clusters, in fact, these points are the same average points of each cluster.  Assigning each sample data to a cluster that data has the smallest distance to the center of that cluster [12]. In the simple form of this method, first, the points are selected randomly as much as needed clusters. Then, the data is assigned to one of these clusters according to the similarity and so, new clusters are obtained [46].  New centers can be calculated for them in each of iterations by repeating the same steps and averaging of data and again the data can be attributed to new clusters [14].
The important steps of this algorithm are summarized as follows [47]: 1) First, k members randomly are selected as the number of clusters among the n members (k is the number of clusters).
2) Z j vector is calculated based on Equation (5) which represents the center of each class C j . It should be noted that relation (6) is used to calculate the center of each class during solving and usually, k samples are randomly selected at the start of the algorithm and are considered as the center of each class [13]. , , , 5) Minimize the objective function of Equation (6) and find the proper classification on the M set with the number k of classes.
And a software has been introduced by the author to speed up the operation above [13] [48].

Concentration-Volume Fractal Method
The C-V fractal model, which was proposed by Afzal et al. (2011) [49] for division of mineralized zones and barren host rocks in porphyry deposits, can be addressed as Equation (7):

Drilling and Sampling
In this study, information of 57 boreholes which are drilled in Sheytoor deposit were used. The depth of the boreholes was between 76 and 475 meters. 1550 samples from cores were systematically collected, the average sampling distance was about 2 meters. All samples were analyzed using ICP-AES method [35]. You can see the geolocations of boreholes in Figure 6.
The statistical parameters of the Fe, P and S elements such as mean, standard deviation, variance and etc. were calculated and are shown in Table 1 [37].

Censored and Outlier Data
The explained method in Section 3. considering these values, which was 24.78% and 9.91%, respectively. According to the graph that is given in Figure 4 and number of samples (1550 samples), the amount of (g) parameter, was estimated "4.4". So using doerffel relationship (8): The highest Fe content in the analysis data is 65.5% which this value is less than outlier limit (68.38%). Therefore, it does not count as an outlier value. So some value is not removed from the raw data as outlier data.

Raw Data Computing
In this section, the frequency distribution chart (histogram) is drawn to examine the various characteristics of our statistical society. The statistical society of the study is a set of geochemical information in the form of iron element concentration derived from systematic sampling of exploratory boreholes which are drilled in the Sheytoor mining area. The histogram of our statistical society is shown in Figure 7. The purpose of this study is geochemical modeling of iron interfering elements (Phosphor and Sulfur) and to identify the behavior of these elements in relation to the Fe.
The statistical parameters of the Iron, Phosphorus and Sulfur elements such as mean, standard deviation, variance and etc. were calculated and are shown in Table 1.

3D Geochemical Model of Fe Element
In this Model, by using the "Inverse Distance Weighting (IDW)" interpolation

Geochemical Behavior of Disturbing Elements
In this section, using the K-Means clustering method, how the elements behave to each other is discussed. The elements associated with the mineralization of iron in Sheytoor area, are sulfur and phosphorus, which can be effective in the processing and extraction of pure metal. Therefore, in this section, the geochemical behavior of iron, phosphorus and sulfur in relation to each other was investigated in pairs.

Geochemical Behavior of Iron and Phosphorus
In order to identify the geochemical behavior of iron and phosphorus elements In order to find the best clustering according to the mentioned profiles, the graph of the mean value of the utility function (S(i)) of the average for the number of classes (K = 3 to K = 10) is given in Figure 10.
According to the above diagram, the highest value of the utility function, which indicates the optimal clustering, occurred in the number of clusters 3. Therefore, we use this clustering to investigate the geochemical behavior of phosphorus with respect to iron. The centers of the designated categories for K = 3 are plotted in Figure 11.
Based on this classification, according to the diagram presented in Figure 11, with the increase of iron concentration in Sheytoor deposit, the amount of phosphorus concentration also increases. The behavior of these two elements in relation to each other is nonlinear. The fitted curve with regression coefficient R 2 = 1 is a quadratic curve whose equation is given in Equation (9).

Geochemical Behavior of Iron and Sulfur
In order to investigate the geochemical behavior of iron element compared to sulfur element in Sheytoor deposit. To select the best classification, the graph of the mean value of the utility function (S(i)) versus the number of classes for iron and sulfur elements and clustering profiles for K = 3 to K = 6 are presented ( Figure 12 and Figure 13). Because the value of the utility function was always decreasing for categories with more than K = 10, only K = 1 to K = 10 was considered.

3D Geochemical Model of P and S Elements
The three-dimensional model of mineralization and dispersion of phosphorus and sulfur elements in exploratory boreholes in order to provide a better view of mineralization in the region is given in the figure below ( Figure 15). Using the Kriging interpolation method in confined space to exploratory boreholes, the estimated titers of phosphorus and sulfur elements were calculated at points where no borehole had been drilled. The mineralization of phosphorus and sulfur elements in the area of Sheytoor deposit was modeled in three dimensions. The 3D geochemical models of P and S elements are presented in Figure  16.

Phosphorus Mineralization
In order to identify the mineralized areas of phosphorus along with the mineralized areas of iron element in Sheytoor deposit, using the three-dimensional model of the dispersion of phosphorus concentration in the area obtained by kriging estimator, and by applying the concentration-volume (C-V) fractal method, fine communities were identified in the area. The threshold of background and anomaly for the element P were also calculated.
As shown in Figure 17, three mineralization communities have been identified for the phosphorus element. The threshold value of each of these communities is presented in Table 2.
The phosphorus mineralization enrichment community is shown in Figure   18. This 3D model is made taking into account the cutoff grade of 1.42%.

Sulfur Mineralization
In order to match the mineralization zones of iron element with high-sulfur and low-sulfur mineralization zones, a three-dimensional model of mineralization and dispersion of sulfur element was constructed using Kriging interpolation  Based on C-V fractal diagram (shown in Figure 19), three communities of Sulfur mineralization are separated. The characteristics of these mineralization communities are presented in Table 3.
The three-dimensional model of mineralization of the anomalous sulfur mineralization community is presented in Figure 20. This model is made taking into account the cutoff grade of 1.07%.

Conclusions
 The region is located in yazd province (central tectonic zone of iran). Statistical and geostatistical studies were carried out based on core data of Sheytoor iron deposit which is located near Bafgh city. The main ore deposit is as apatite-magnetite mineralization.  In order to identify and replace the outlier data, the doerffel method was used. It was found that there is no outlier data among the geochemical data of the iron element obtained from the analysis.  Using Kriging and Inverse Distance Weighting (IDW) interpolation methods, 3D block models were created and presented.  Figure 19. Separation of sulfur mineralization grade communities using C-V fractal method.  Open Journal of Geology  Using the K-means clustering method, the concentration data of phosphorus and sulfur elements in relation to the element iron were classified. Then the optimal clustering was determined and based on the concentration of elements in the centers of the classes, the geochemical behavior of phosphorus and sulfur elements towards the iron element were modeled.  After mineralization and dispersion of phosphorus and sulfur concentrations were modeled in three dimensions, mineralization communities were separated for each element using C-V fractal method.  The enrichment zone threshold for phosphorus and sulfur are 1.42% and 1.07%, respectively. These communities were modeled three-dimensionally to adapt to iron mineralization in the Sheytoor iron ore mine.  Mining planning can be reconsidered using three-dimensional models of mineralization of phosphorus and sulfur elements, which are considered as deleterious elements in the process of iron concentrate production and, of course, later stages such as steelmaking. It is suggested that mineral extraction specialists continue geometallurgical studies based on the information obtained from this study.

Conflicts of Interest
The authors declare no conflicts of interest regarding the publication of this paper.