Evaluation of heavy metal contamination using cokriging geostatistical method (case study of Abteymour oilfield in southern Iran)

One of the most significant sources of soil contamination is the industries involved in the discovery, extraction, and use of petroleum and gas resources. One of the largest oil fields in Iran is the Abteymour oilfield, which is situated in the agricultural areas and Karun river flood plain. This study was carried out to analyze and monitor soil contamination by heavy metals and to create a map of the spread of contamination in the vicinity of the Abteymour oilfield, taking into account the significance of soil pollution in such a region. Thirty-three samples from the local surface soils were used in this study, and in addition to testing the ICP-MASS device's ability to detect the presence of heavy metals and other key elements, some of the soil's physical and chemical characteristics were also determined. After drawing the variograms and determining the appropriate fitting model, the cokriging method was used to study the spatial distribution of Al, Ti, Sr, Cr, Ni, V, Cu, Li, and Pb heavy metals and Na and S elements. The results of the heavy metals distribution in the studied area showed that the distribution of the studied elements is affected by geological conditions and human activities, including the petroleum industry and agriculture activities.


Introduction
The presence of heavy metals in soils can be brought on by anthropogenic activities or natural processes like weathering. Heavy metals in the soil may accumulate more due to human activity (Yalcin et al. 2007). Various causes, including sewage sludge, chemical fertilizers, wastewater leaks and infiltration from various industries and mines, petroleum facilities, landfills for various sorts of waste, fossil fuels, landfilling of different forms of waste, and so on, cause heavy metal contamination of soil and water resources (Wei and Yang 2010;Mahmoudi et al. 2013). One of the most significant sources of soil contamination is the industries related to the discovery, extraction, and use of oil and gas resources. Due to their vital role in delivering energy and creating raw materials for many other industries, the oil and gas industries are particularly significant. These sectors are among the most important because of this.
The heavy metal contamination in the processes of drilling and extraction of petroleum resources are originated from drilling fluid and drilled formations, as well as the extracted petroleum and gas. During the drilling operation, the drilling fluid and the cuttings of various drilled formations, which are also impregnated with petroleum in some cases, are accumulated in the pits around the drilling site. If environmental considerations are not considered in the accumulation of these wastes, over time, various contaminants will penetrate into the soil environment and then contaminate groundwater resources over time. As well as, the wind and surface water flow spread of contamination in a wide region. The activity of refineries causes further heavy metal contamination via soil contaminated with petroleum products. Furthermore, soil contamination with petroleum products is workable during petroleum transfer processes. Along with this petroleum contaminates, heavy metals such as nickel (Ni), vanadium (V), lead (Pb), iron (Fe), cadmium (Cd), zinc (Zn), and copper (Cu) penetrate the soil as well (Adesina and Adelasoye 2014 Several studies have been conducted in relation to soil contamination with heavy metals; Instances of which are listed in the following: The environmental effects of the drilling mud system on the desert region of Algeria were examined by Ghazi et al. (2011). This study looked at the pollution of heavy metals caused by various drilling waste management techniques, including barium (Ba), antimony (Sb), arsenic (As), aluminum (Al), zinc, and aromatic hydrocarbons. To carry out this study, the researchers employed the Life Cycle Assessment methodology. The findings revealed that the majority of the heavy metal concentrations were higher than the recommended levels, and that stabilizing/solidifying the waste already in place online was the most effective way to reduce contamination in this area. Karbassi et al. (2015) measured the concentration of eight heavy metals including Cu, Ni, V, Cd, Pb, Zn, cobalt (Co), and molybdenum (Mo) in the soils around the oil fields of Ahvaz and compared them with the average of the Earth's crust. The researchers also used two indicators of geo-accumulation (Igeo) and enrichment factor (EF) to check the level of pollution. The results showed that all the studied metals had higher values than the average crust. Also, based on the Igeo, it was found that the intensity of pollution in the area changes from unpolluted to highly polluted. The researchers suggested using phytoremediation method to remove these elements. Doyi et al. (2018) studied the distribution of heavy metals, including manganese (Mn), Cd, Ni, Pb and As in the Tano Basin in Ghana using Igeo, risk assessment models and multivariate data analysis techniques. Based on the Igeo, it was found that moderate to strong Pb accumulation levels while all the other metals show uncontaminated to moderate levels. The Igeo showed that the accumulation of lead was moderate to strong, while other metals show uncontaminated. Additionally, human activities were introduced as the source of Cd and Pb pollution in this area. Miao et al. (2019), was studied the distribution of As, Cd, Zn, Pb, Cu, chromium (Cr), and mercury (Hg) in the soils of the Yellow River Delta in China using environmental indicators and spatial distribution. The results showed a nearly toxic concentration of mercury and further established that the industrial and petroleum activities had caused more contamination than agriculture in the region. Qaiser et al. (2019) investigated the concentration of heavy metals (Ba, Pb, Cr, Cd, Zn, Mn, and Ni) in soil and drilling mud waste in a region of Pakistan using environmental indicators. The results showed that Pb and Ba concentrations were higher than the standard and had caused pollution. Therefore, statistical results indicated the same origin of heavy metals in this region.
There are several methods to investigate and monitor heavy metal contamination in the soil. Accordingly, geostatistical methods are efficient tools for modeling the distribution of contaminations and determining the potential origins over wide regions. Due to the spatial limitations and sampling problems to determine polluted areas and those exposed to contamination, the use of geostatistical methods has come to the attention of many researchers and experts in recent decades (Lu et al. 2012;Lv et al. 2015;Mihailovic et al. 2015;Zhou et al. 2016;Liang et al. 2017;Hung et al. 2019;Chuna et al. 2019).
For estimating variables with temporal and spatial variation, there are numerous geostatistical techniques. The estimation of the weight factor, which is generalized to the observed points around the desired point, is where these methods diverge most from one another. Among the geostatistical methods, the Cokriging method has a high potential to estimate the spatial distribution of heavy metals in the soil. Furthermore, it has been suggested as a suitable method for interpolation and preparation of contamination maps (Burgess and Webster 1980;Juang et al. 2001;Khosravi et al. 2018;Wang et al. 2018;Cunha et al. 2019).
In light of the significance of soil contamination in industrial areas, particularly oil fields, this study was carried out to look into and track soil contamination by heavy metals and create a map of contamination distribution in the vicinity of the Abteymour oilfield in southwest Iran. This oilfield is located in an agricultural region. Therefore, it is very important to identify the type and amount of possible contamination to prevent it from entering the food cycle. In addition, the Karun River, which is the longest and highest water flow river in Iran, passes through the middle of this oilfield. Soil contamination in the area of the oilfield will lead to river contamination through surface water flow. Considering that the Karun River is fed from underground water sources, the penetration of contamination into these sources can also cause river contamination. This river finally leads to the Persian Gulf. In this research, geostatistical cokriging method will be used to model the spatial distribution of heavy metals contamination and investigate the effect of petroleum activities on the distribution of contamination in the Abteymour oilfield. Just as there are multivariate methods for estimation in classical statistics, in geostatistics, interpolation can be done with more accuracy using Cokriging methods that work based on correlation between data (Hassani Pak 2013). In some instances, the data set includes one or more variables that spatially depend on the primary variable and reveal important details about it (Deutsch and Lewis 1992). The use of cokriging method has grown in soil and environmental studies due to factors like its capacity to use the data of secondary variables to improve estimation accuracy, mutual correlation between measured variables, economic, temporal, and technical issues with measuring some variables, increasing processing power, and the existence of diverse and readily available software (Alvares et al. 2013).

Geological setting
Abteymour oilfield, 23 km in length and 6 km in width, is located 68 km southwest of Iran (Fig. 1). The oilfield is located in the Zagros zone and the northern Dezful embayment. The structural-sedimentary Zagros zone is one of the biggest morphotectonic belts of Iran, from the northeast to the southwest, include the thrust zone, the simple folded zone and the Dezful embayment. The Dezful embayment includes most of Iran's oilfields. Structures with NW-SE trend are the characteristics of this basin. Structures with N-W trend such as Khorramshahr and Darkhoein oilfields are also seen in this basin. Figure 2 displays the area of the Abteymour oilfield with the NW-SE trend in relation to the geological formations, which are composed solely of Quaternary alluvial deposits.

Sampling method
Due to the limited access to the whole area of Abteymour oilfield for sampling, in this research, in order to better cover the area, a systematic random sampling method was used. As well as, according to the limitation of the number of point-pairs in the variogram drawing, which should not be less than 15 point-pairs (Hassni Pak 2013), collection at least 33 samples (considering possible errors) was considered. Therefore, the area of Abteymour oilfield was gridded, and the location of 33 samples (one sample from each cell) was determined randomly (Fig. 3). Samples were collected from surface soils (0-20 cm depths). The intended parameters in this study are measured according to Table 1.

Methodology
Geostatistical methods are helpful tool in studying variables that have a spatial structure (Wackernagel 2003;Hassni Pak 2013); therefore, at first, the spatial structure of the data must be determined, and the data analysis should be down after confirming the structure. Spatial structure can be identified and investigated by variogram. The variogram shows differences between the sampled values as a function of distance between point-pairs. Furthermore, concepts of continuity and homogeneity of samples are also examined by variogram (Cressie and Hawkins 1980). Based on the idea that the similarity of a phenomenon properties decreases with increasing distance, the variogram measures the continuity degree (dependence or autocorrelation) between points (Chrisman 2002). Spatial autocorrelation depends on the distance between point-pairs and changes with the variation of distance. Stable spatial variation of continuity is called isotropy. In addition to distance, continuity may also change with direction. The effect of a direction change on a variogram is described as anisotropy (Chilès and Delfiner 2012). Anisotropy is usually not a deterministic process that could be described simply by a mathematical formula. In fact, anisotropy is a random process that causes more continuity in a specific direction. These variations can be the result of some unknown and unmeasurable physical processes that are modeled as a random process with directional continuity. In addition to anisotropy of the variogram, general trends are also effective in the modeling process. General trends are an essential process that influences all dimensions of a modeling process. Trends in environmental data show that the value of a variable is not a constant in different directions and varies based on the distance of the points to each other. The trend can be identified by a three-dimensional diagram and regression equations. If the line fitting on the data distribution in the three-dimensional diagram is flat and fixed, it signifies the absence of a trend. If the trend deters accurate modeling of the estimated surface, it should be removed by appropriate regression equations.
Nugget effect (C 0 ), sill (σ 2 ), and range (R) are aspects of a variogram that are crucial in geostatistical analysis. The nugget effect is represented by the variogram's value at the origin. Given that the spatial organization can be more accurately modeled, this number is nearer to zero. Sill is represented by the variogram curve's upper bound. The maximum distance at which the examined variable has spatial correlation is displayed by the range, also known as the range of influence or the range of correlation. It goes without saying that a greater range denotes a wider spatial organization. While anisotropic variables have a range with two directions, the maximum and minimum, isotropic variables have an identical range of impact in all directions (Rouhani et al. 1996). A useful measure of association is the nugget to sill ratio (correlation ratio), which is reported as a percentage. If this value is less than 25%, the variable has a high spatial correlation and the role of the structural component is greater. If the value is between 25 and 75%, the correlation is moderate. Otherwise, the correlation of location is weak. In such cases, the use of geostatistics will not be very useful (Cambardella et al. 1994;Vierira and Paz-Gonzalez 2003).
The most appropriate theoretical models must be fitted to the experimental variogram before its application in modeling. As the number of probability distribution models in classical statistics is limited, a small number of well-known models are commonly used for variograms, which include (a) bounded models (such as pure nugget effect, spherical, exponential, and Gaussian), and (b) unbounded models (such as linear, Davison, and parabolic) (David 1977;Journel and Huijbregts 1978).
Variogram-based modeling employs a variety of techniques, mainly falling into the deterministic and stochastic categories (Dubrule 2003). An average, homogeneous, and smooth solution is estimated using deterministic approaches, one of which is the kriging method. Contrarily, probabilistic approaches mimic a large number of likely but diverse and chaotic events; one example of such a method is Gaussian Sequential Simulation (GSS).
Kriging is one of the most extensively used and comprehensive estimation methods, in which most important conditions relevant for estimation are reflected. In this method, which is based on the weighted moving average, the variances of the prediction error are minimized. Therefore, this method is introduced as the best punctual or block unbiased linear estimator (Oliver and Webster 2015). In addition to the estimated quantity, this estimation method calculates the error associated with each estimate; therefore, the error distribution is also determined. There are numerous types of kriging; the selection of each of them for different purposes depends on the conditions of the studied variable (including having a normal distribution and having a trend or lack it). The normality of variable distribution is one of the main assumptions of kriging, although they are rarely investigated in case studies (Hengl 2007). Kriging models are also divided into linear and nonlinear categories (Sarma 2001). In this division, normality is a necessary condition for using linear kriging methods (simple and ordinary kriging) (Hengl 2007). Ordinary kriging is robust even when it deviates from basic assumptions; for this reason, it is the most popular method. If the data distribution is non-normal, nonlinear kriging methods must be used or the data converted  to normal via conventional methods (Rouhani et al. 1996;Webster and Oliver 2008;Sarma 2001). If the trend is not removed, the universal kriging method should be used for modeling. Utilizing a variable's linear combination with other variables is the most popular and straightforward way to estimate it in a particular location (Alijani et al. 2013). Cokriging is an extension of kriging in which the interpolation accuracy of the examined variable is enhanced by using extra observable variables (known as covariates, which are frequently connected with the studied variable). Co-located samples, additional target sites, or both can be used to measure the secondary variable. In reality, by utilizing both the spatial organization of the data and the interactions between several variables, the cokriging approach extends the kriging method to multivariate interpolation. In this study, the studied factors were modeled using the conventional cokriging approach.
Output evaluations are necessary to confirm the accuracy of each interpolation method. With several methods for this purpose, cross-validation is regarded as the most significant method. Based on this method, measured and estimated values of a variable are compared through different indices. Root-mean-square error (RMSE) and mean error (ME) are the most significant of the different criteria for this end and have been used in this study as well. These criteria are determined according to the following equations, and values closer to zero show lower errors (Webster and Oliver 2008).
Z (x i ): measured values of the variable in x i , Ẑ x i : estimated measured values of the variable in x i , N: number of data.

Results and discussion
In order to evaluate the distribution of heavy metals in the Abteymour oilfield, metals that have values higher than the background concentration or are important in terms of regional geology were investigated. Accordingly, in addition to the Al, Cr, Ni, Cu, Pb, V, strontium (Sr), titanium (Ti), and lithium (Li) metals, these elements contain sulfur (S) and sodium (Na) elements. Table 2 shows the descriptive statistical indices for the studied elements and used variables in the cokriging method. Based on the results of the Kolmogorov-Smirnov (K-S) test to determine the normality of data, studied parameters, except Na and S, are normal (Table 3). These two parameters were converted to normal data by a logarithm. Table 4 displays the correlation of the studied parameters. The second step in modeling is to identify the general trends of the data. According to Fig. 4., the Na, Al, Sr, Ni, V, Cu, Li, and Pb elements have trends in NE-SW and NW-SE directions, and the Cr, S, and Ti elements have the N-S and E-W trends. According to the resulting curves, the quadratic equation was used to remove the trend of the variables. Accordingly, the conventional cokriging method was used for geostatistical analysis due to the normality of the studied parameters and the elimination of the trend. Due to the normality of the studied elements and removal of the general trend, the ordinary cokriging method was used for geostatistical analysis. Figure 5 and Table 5 show the result variograms for the studied elements using the cokriging method. Based on the correlation in the table, all the elements have a spatial structure. Similarly seen in Table 5, the Al, S, Ti, V, and Li elements have anisotropy, while the Na, Sr, Cr, Ni, Cu, and Pb elements are isotropic. Figure 6 shows the distribution of the studied elements based on cokriging modeling. As explained earlier, all elements have trend(s), which can be clearly seen in the map of estimated surface. This general trend could be the result of changes in natural or artificial factors. Agricultural and industrial activities are two of the important artificial factors.
1. Al According to the map, Al concentration is the highest in the S-SW and NE direction. Due to the positive correlation of Al with V and Ni as well as the similarity of their distribution, it can be concluded that these elements have a common origin. However, the interference of different agricultural and industrial activities in the studied area complicates the conclusion regarding the origin of these elements. 2. Na The concentration of Na element is high throughout the region, which, in addition to the environmental and geological conditions of the region, can be caused by the increase in agricultural (due to excessive use of chemical fertilizers) and industrial activities. The Karun River, long used for irrigation purposes in the region, has undergone an increase in salinity due to reduced discharge and increased inflow of various wastewater. Consequently, the river's salinity could cause soil salinity. 3. Ti The concentration of Ti is highest in the NE and SW parts of region. This element has a positive correlation with Ni and V elements, which are crucial elements in petroleum and gas industry. This correlation can indicate their common origin.

S Concentration of S is highest in the central and SE
parts of region. Agricultural activities (use of chemical fertilizers) and waste from drilling oil wells are among the most important factors affecting the concentration of this element. Total S in Iranian heavy crude oil amounts to 2.28 percent, while it is 1.58 percent in Iranian light crude oil. 5. Cr Concentration of Cr has increased significantly in the E and SE parts, as well as a small part of the NE. Sewages from various petroleum and agricultural Due to the use of chemical and animal fertilizers, agricultural activities can increase this element in the soil. 7. V The Ni, Al, Li, and V elements have a similar trend.
These elements have the highest values in the NE, E, and SW areas, which could denote their common origin. Altogether, mentioned elements are the components of crude oil. 8. Ni The highest Ni concentrations are observed in the NE, W, and SW parts of the region. Human activities such as drilling oil wells and agriculture are the most influential factors in the rise of Ni concentrations in the region. Accordingly, n regions with shale or marl shale rock, the Ni concentration is more influenced by the lithology. Based on the specified limits for the  concentration of Ni in Iranian soils (50 mg/kg), 85% of the samples have a concentration that exceeds the permissible limit. 9. Li The highest Li concentration is seen in the NE, SW, and E parts of region due to the presence of this element in the chemical structure of urea, phosphate, and potash fertilizers, excessive use of chemical fertilizers in agricultural lands has increased its concentration in the soil as well as the urban and rural wastewater have increased Li concentration in soil.

Cu
The highest Cu concentrations are observed in the central of the region. The overload of oil drilling, excessive use of chemical fertilizers in agricultural lands and sewage are the main reasons for the increase in Cu concentration. Due to the presence of this element in the chemical structure of urea, potash, and phosphate fertilizers, improvident use of chemical fertilizers in agricultural lands has increased its concentration in the soil. Disposal of household wastewater, household waste, and wastewater from drilling oil wells are other factors that increase the Cu contamination in the Abteymour oilfield. 11. Pb Highest Pb concentration is observed in the N and SE parts, as well as a small part of the SW region. Pb has diverse applications and penetrates water and soil sources from different origins. This element is found in the chemical structure of urea, phosphate and potash fertilizers, domestic wastewater, and vehicles smoke.
Drilling, extraction, and exploitation processes used in the oil and gas sector are effective contributors to the contamination of an environment with heavy metals. The release of drilling fluid, sewage, and trash, as well as oil leakage, all have an impact on the detrimental effects of these activities (Asia et al. 2007). Environmental dangers are caused by heavy metal contamination, which persists in the environment for many years (Ebenezer and Eremasi 2012).
In recent years, there have been many concerns regarding the effect of drilling fluids and their additives in to drilling fluid to increase the fluid weight to control well pressure. Heavy elements affiliated with barite that have environmental consequences include Pb, Zn, Cu, Hg (Crecelius et al. 2007), and Cd (Nelson et al. 1984), which are insoluble as sulfide salts (Leuterman et al. 1997). Due to low bioactivity, these elements are non-toxic (Schaanning et al. 2002) and are released only under highly acidic conditions. Drilling wastes includes cuttings that may also be contaminated with crude oil, and drilling fluids are usually discharged into pits near drilled wells. Due to these pits are opened, the contamination could spread to the surrounding environment under climatic conditions. Accordingly, Cr, Zn, Sn, Ni, Cu, Cr, Co, Ba, and Sr in drilling mud wastes are determined to be lower than the permissible limits while Hg and As value.
Abteymour oilfield is located among agricultural lands (Fig. 7), and the drilling and extraction activities of petroleum ( Fig. 8) have caused discontent among farmers. Failure to comply with environmental standards in the drilling and storage of industrial materials and petroleum wastes that are deposited in inadequate pits and the absorption of these toxic substances in the soil and groundwater has severely damaged the region. Cr, Ni, V, and Mg elements are mostly of natural origin and components of crude oil. These elements perhaps penetrated the drilling fluid and then the soil of region as a result of the drilling operations. Metals such as As is used for disinfection and Ba is used as the main component of barite additive. Cd and Cr elements are included in the composition of drilling fluid in the form of chromium lignosulfonate additives or impurity or other components of drilling fluid. In addition to the drilling fluid composition, high levels of Pb are due to the equipment, drilling rigs, and fuels used. Al is the fundamental component of silicate minerals, an element that is released into the sediments of the Earth's surface as a result of weathering. The distribution of Al is mainly controlled by the sediment gradation, particularly of fine particles such as silt and clay (Zhou et al. 2014).

Conclusion
The most important results of this study are as follows: -Investigation of heavy metals in the Abteymour area showed that Li, Pb, Na, Al, Sr, Ni, V, Cu elements have a trend in NE-SW and NW-SE directions, while Cr, Ti and S elements have a N-S and EW trends. Quadratic equations were used to remove the trend. -The resulting variograms that the studied elements have a spatial structure. The Na, Sr, Cr, Ni, Cu, Pb elements have the smallest range and isotropic distribution, while the Al, S, Ti, V, Li elements have the largest range and anisotropic distribution. -The Al distribution reflects the improvident use of chemical fertilizers in agriculture, herbicides and pesticides, as well as the wastewater from fish ponds. -The Na element has a high concentration in the region.
The geological structure and agricultural activities (due to the presence of Na in the chemical structure of urea, phosphate and potash fertilizers) have increased the concentration of Na in the soil. -The Ti distribution is influenced by the by petroleum activities in the study area. -The S element is influenced by both natural and human factors in the central and southeastern prats of region. Agricultural and industrial activities (including drilling for oil wells), cause S to penetrate the soil and the environment. -The Cr distribution is the highest in the eastern, southeastern, and a small portion of the southwestern areas. The parent rock and agricultural activities has increased this element in the region. -The distribution of Sr element is affected by human and agricultural activities (due to excessive use of chemical and animal fertilizers). -The V, Ni, Ai, and Li elements have similar distribution trends that are influenced by their common origin. These elements are the components of crude oil. -Due to the Li element is found in the chemical structure of urea, phosphate, and potash, improvident use of chemical fertilizers in agricultural lands has increased Li concentration in the soil. -The highest Cu concentration is observed in the central parts of the region, where oil wells are concentrated. The discharge of domestic wastewater, urban waste, as well as wastewater and waste from drilling oil wells are among the factors that increase Cu contamination in the Abteymour oilfield. -The highest concentration of Pb element is observed in the north, southeast, and a small part of the southwest of region. Fish ponds, roads, and urban contaminations are widespread in these areas.