Environ Sci Pollut Res (2016) 23:2758–2769 DOI 10.1007/s11356-015-5507-2
RESEARCH ARTICLE
Comparison of interpolation methods for the estimation of groundwater contamination in Andimeshk-Shush Plain, Southwest of Iran Rouhollah Mirzaei 1 & Mohamad Sakizadeh 2
Received: 7 August 2015 / Accepted: 27 September 2015 / Published online: 7 October 2015 # Springer-Verlag Berlin Heidelberg 2015
Abstract Selection of appropriate interpolation methods for the conversion of discrete samples into continuous maps is a controversial issue in the environmental researches. The main objective of this study was to analyze the suitability of three interpolation methods for the discrimination of groundwater with respect to the water quality index (WQI). The groundwater quality data consisted of 17 variables associated with 65 wells located in Andimeshk-Shush Plain. Three spatial interpolation methods including ordinary kriging (OK), empirical Bayesian kriging (EBK), and inverse distance weighting (IDW) were utilized for modeling the groundwater contamination. In addition, different cross-validation indicators were applied to assess the performance of different interpolation methods. The results showed that the performance differed slightly among different methods, although the best performed interpolation method in this study was the empirical Bayesian kriging. Among the interpolation methods, IDW with weighting power of 4 estimated the most contaminated area, while OK estimated the lowest contaminated area. The weighting power of IDW had a significant influence on the estimation, meaning that the estimated contaminated area was increased when a greater weighting power was selected. The subtraction results indicated that there are slightly spatial differences among the contamination assessment results. Results of both standard deviation (SD) and coefficient of variation Responsible editor: Marcus Schulz * Rouhollah Mirzaei [email protected] 1
Department of Environment, Faculty of Natural Resources and Earth Sciences, University of Kashan, Kashan, Iran
2
Department of Environmental Sciences, Faculty of Sciences, Shahid Rajaee Teacher Training University, Tehran, Iran
(CV) also showed that uncertainty was highest in the southern part of the study area, where the distribution of wells were more intensive than that of the northern part. Keywords Groundwater contamination . Uncertainty . Inverse distance weighting . Empirical Bayesian kriging . Iran
Introduction Groundwater is a precious water resource which is used in different fields such as drinking, irrigation, and industrial processes (Nampak et al. 2014; Neshat et al. 2015). As a convenient substitute, groundwater also provides additional water supply to compensate for insufficient surface water resources. As such, groundwater should be protected from contamination by implementing effective water resource planning and management (Neshat et al. 2015). Iran is one of 27 countries that are likely to face increasing water shortage crisis by 2025 unless action is taken to reduce current water consumption (Nabi Bidhendi et al. 2007). Most of the groundwater resources in Iran are used as drinking water besides agricultural and industrial purposes. So, the quality of groundwater is of great importance. Low quality water supply can cause health problems; therefore, determining the quality of water is an important issue. The estimation of spatial distribution of contaminated groundwater is very important in the health risk assessment (Lee et al. 2007). However, only a small proportion of in situ data can be analyzed in a field investigation, due to time and cost constraints. Consequently, sparse measured data contain considerable uncertainty (Liu et al. 2004). Therefore, mapping the spatial distribution of groundwater contamination requires spatial interpolation methods. Surface interpolation methods in a geographic information system (GIS) are very powerful
Environ Sci Pollut Res (2016) 23:2758–2769
tools for predicting surface values (Gong et al. 2014). The discipline of geostatistics provides very useful techniques for handling spatially distributed data such as groundwater contamination (Arslan 2012; Nas and Berktay 2010). By identifying spatial patterns and interpolating values at unsampled locations, interpolation analysis can play a vital role in the sustainable management of groundwater systems by providing estimated input parameters at regular grid points from measurements taken at random locations (Arslan 2012; Kumar 2007). Many authors have emphasized the role of geostatistics in the management and sustainability of regional water resources (Arslan 2012). The most commonly used interpolation methods are inverse distance weighted, kriging, and cokriging interpolations (Gong et al. 2014). However, it has been inconsistent concerning a potential superiority of one interpolation method over another (Gong et al. 2014). Consequently, interpolation methods such as inverse distance weighting (IDW) and kriging have been extensively used in groundwater investigations and contamination mapping (Arslan 2012; Gong et al. 2014). Interpolation accuracy is related to the precise definition of the contaminated area and its boundaries. Consequently, this directly affects the accuracy of contamination assessment. There are a lot of studies on the performance of the spatial interpolation methods for different objectives, but the results are not conclusive. Through interpolation of different data set, some of the authors found that the kriging method has the lowest prediction bias among other methods (Joseph et al. 2013; Liu et al. 2014; Yasrebi et al. 2009), while some other authors showed that IDW was the best method (Gong et al. 2014; Keblouti et al. 2012). Moreover, there are other studies in which different interpolation methods had outperformed such as TPS (Gumiere et al. 2014), ANN kriging (Dai et al. 2014), and Bayesian kriging (Plouffe et al. 2015). Interpolation methods all have a smoothing effect, which underestimates the local high values and overestimates the local low values (Journel et al. 2000; Xie et al. 2011). This smoothing effect leads to bias (e.g., underestimation or overestimation) in contamination assessment and has an effect on relevant environmental decision-making (Goovaerts 2000; Xie et al. 2011). It is essential to minimize this bias in a contamination assessment caused by interpolation methods and to understand the uncertainty of groundwater contamination assessment introduced by interpolation error and differences in contamination assessment between various interpolation methods (Xie et al. 2011). Considering the importance of estimating groundwater contamination, this study was undertaken to compare the performance of three interpolation methods, including IDW, ordinary kriging (OK), and emperical Bayesian kriging (EBK), in the estimation of groundwater quality in Andimeshk-Shush Plain using the leaveone-out cross-validation.
2759
Material and methods Study area Andimeshk-Shush Plain covers an area of approximately 1100 km2 in Khuzestan, Iran (Fig. 1). This plain is located in the north of Khuzestan (32° 00′–32° 35′ N; 48° 10′–48° 25′ E). The local economy depends largely upon farming. Tourism and manufacturing also contribute to the area’s economy. Andimeshk-Shush Plain is comprised of a succession of Dez and Balarood seasonal river deposits interspersed with minor silt in Andimeshk Plain in addition to sand and gravel deposits of Dez and Kharkhe Rivers interspersed with major clayey silt lenses in Shush Plain. Various lithological units range from Cenozoic (Pliocene) to Quaternary Periodare are located within the Andimeshk-Shush Plain. Quaternary-age deposits consist of alluvium, which have been made up of loose, interlayer clay, silt, sand, and gravel. The thickness of the alluvium is about 200–300 m. Andimeshk-Shush aquifer is the primary source of groundwater, supplying approximately 100 % of the total drinking water for about 180,000 people residing in this area. Annual precipitation is approximately 270 mm from which more than 80 % occurs during the December–April period. In comparison, annual potential evaporation is about 1670 mm that is six times higher than that of the annual precipitation. Farms occupy over 70 % of the study area, in which the main agricultural crops are wheat, corn, and sugarcane. Nearly more than 75 % of farmlands in the study area are irrigated by the surface water from Dez irrigation network and the rest (less than 25 %) are irrigated with groundwater. Groundwater quality parameters and water quality index calculation The groundwater quality data used in this study included 17 variables (EC, TDS, turbidity, pH, total hardness, Ca, Mg, sulfate, phosphate, nitrate, nitrite, fluoride, chloride, Fe, Mn, Cu, Cr(VI)) which have been sampled periodically by Andimeshk Health Network and Iran Ministry of Energy from 65 wells for 8 years between 2006 and 2013. Minimum, maximum, average, and standard deviation of these parameters have been given in Table 1. To calculate the water quality index (WQI), the method proposed by Horton (1965) and followed by many researchers Patel and Desai (2006) was used. For computing WQI, three steps were followed as explained by Rupal et al. (2012). We used annual mean of each parameter for WQI calculations. In the first step, a weight (wi) was assigned to each parameter according to its relative importance in the overall quality of water for drinking purposes. The highest weights were assigned to parameters with the most detrimental effect on the health of groundwater consumers, while the parameters with the least issue in
2760
Environ Sci Pollut Res (2016) 23:2758–2769
Fig. 1 A sketch map of the study includes the distribution of sampling wells in Andimeshk-Shush Plain
this regards took the lowest weights. In the second step, the relative weights (Wi) were computed according to this equation: Xn w ð1Þ W i ¼ wi = i¼1 i where Wi is the relative weight, wi is the weight of each parameter, and n is the number of parameters. In the third step, a quality rating scale (qi) for each parameter was assigned by dividing its concentration in each water sample by its respective standard and the result multiplied by 100: qi ¼ C i =S i 100
ð2Þ
where Ci is the concentration of each chemical parameter in each water sample, and Si is the associated drinking water standard. The final WQI was determined by the product of Wi and qi as the following: WQI ¼ W i qi
ð3Þ
Inverse distance weighting Inverse distance weighting is based on the premise that the predictions are a linear combination of available data. In this method, the interpolating function is as follows: Xn w i zi Z ðxÞ ¼ Xi¼1 ð4Þ n w i i¼1 In which wi¼di−u where Z(x) is the predicted value at an interpolated point, whereas Zi is the amount at a known point. n is the total number of known points used in interpolation, di is the distance between point i and the prediction point, wi is the weight assigned to point i. Higher weighting values are assigned to those points which are closer to the interpolated point. As the distance increases, the weight decreases and u is the weighting power that impose the amount of weight decrease with respect to the increase in distance (Xie et al. 2011).
Interpolation methods Ordinary kriging In this study, the most widely used interpolation methods, inverse distance weighting, ordinary kriging, and empirical Bayesian kriging were evaluated.
Kriging is a linear estimator meaning that the estimate of the unknown value is a linear combination of the known data
Environ Sci Pollut Res (2016) 23:2758–2769
2761
Table 1 Descriptive statistics of groundwater quality variables in Andimeshk-Shush Plain along with the associated WQIs Groundwater quality variables
Mean (mg/l)
Min (mg/l)
Max (mg/l)
SD
pH
7.91 0.46 37.69 18.70 240.16 77.06 596.88
6 0.02 0 0 0 8 8
8.13 0.95 102.32 85 1125 213 36,545
2.23 0.189 0.189 17.48 203.56 34.37 4528.60
858.69 435.74
121 60
2367 1185
469.21 229.91
F− Cl− NO3− SO42− Ca2+ Mg2+ EC TDS Turbidity
1.77
0
21
3.40
Total hardness Phosphate
340.68 0.16
73 0.034
1200 0.47
172.53 0.073
Nitrite Fe Mn
0.07 0.09 0.2
0.003 0.01 0
0.46 0.27 2.66
0.107 0.077 0.410
Cu Cr(VI) WQI
0.42 0.5 76.61
0.015 0.004 55
2.90 0.076 90.96
0.646 0.015 8.52
values (Xie et al. 2011). The aim of kriging is to estimate the value of a random function, z, at one or more unsampled points or over larger blocks, from more or less sparse sample data on a given support, say z(x1), z(x2),…z(xN), at x1, x2,…xN. This can be shown by Xn ð5Þ w Z xj z * ð x0 Þ ¼ i¼1 i where wj are the weights assigned to the known value of z(xj) and z*(x0) is the estimated value. To ensure that the estimate is unbiased, weights are made to sum to 1 (Xie et al. 2011): Xn w ¼1 ð6Þ i¼1 i
Emperical Bayesian kriging Empirical Bayesian kriging (EBK) is a geostatistical interpolation method that automates the most difficult aspects of building a valid kriging model. In addition to accounting for the uncertainty in the underlying semivariogram parameters, the other main redeeming feature of EBK is that despite the common applied Geostatistical Analyst in ArcGIS10.2, the parameters in the new developed EBK are automatically optimized through a subsetting and simulation process which is implemented by estimating a lot of semivariogram models
instead of a single semivariogram. The prediction in unknown locations in common kriging methods is done through calculation of semivariogram with respect to the known data locations resulting in the underestimation of the standard error of the prediction due to overlooking the uncertainty of semivariogram. On the contrary, EBK uses an intrinsic random function as the kriging model despite the other kriging methods. The other main difference of EBK with that of the other kriging models is that EBK does not assume a tendency toward an overall mean; thus, there is the same chance for large deviations to get larger or get smaller. The following process is followed in EBK. (1) Using the available data, a semivariogram model is estimated. (2) Given this semivariogram, a new value is simulated at each of the input data location. (3) With respect to the simulated data, a new semivariogram model is estimated accordingly. The calculation of a weight for the latest semivariogram according to Bayes’ rule is the next step in this field. The semivariogram estimated in step 1 is used to simulate a new set of values at the input location d u r i n g t h e r e pe t i t i o n o f s t e p s 2 an d 3 . A ne w semivariogram model and its weight are produced given the simulated data. During this process, the predictions and their respective standard errors are produced at the unsampled locations. This process finally creates a spectrum of semivariograms (Krivoruchko 2012a, b). There are two base distributions available in EBK with respect to the utilized multiplicative skewing normal score in this method: empirical and log empirical. Since log empirical only accept positive data values and their predictions are also positive, so it is a good option for water quality indices that are just positive scores. Due to the application of log transformation on our water quality indices, exponential model was applied to interpolate WQI in this study. Comparison of the accuracy of the interpolation methods Cross-validation and validation with an independent data set are the commonly used methods for comparing the interpolation methods. Because the sample size was limited, a leave-one-out cross-validation (LOOCV) was used to estimate water quality index of each well. With LOOCV, one data point (e.g., water quality index in one well) is removed from the dataset, and its value is estimated by the remaining known values. This process is repeated until water quality index values from all the wells are estimated (Xie et al. 2011). In order to assess the performance of different interpolation methods the mean error (ME), mean absolute error (MAE), mean relative errors (MRE), mean squared error (MSE), root mean squared errors (RMSE), Nash–Sutcliffe efficiency
2762
Environ Sci Pollut Res (2016) 23:2758–2769
(NSE), and percentage bias (PBIAS) of the interpolated values were calculated with the following equations: Xn Z * ðxi Þ−Z ðxi Þ i¼1 ð7Þ ME ¼ n Xn Z * ðxi Þ−Z ðxi Þ i¼1 MAE ¼ ð8Þ n 1 X n Z * ðxi Þ−Z ðxi Þ ð9Þ MRE ¼ i¼1 n Z ðx Þ X
MSE ¼
i
Z ðxi Þ−zðxi Þ *
2
n rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X n 2 1 * RMSE ¼ Z ð x Þ−Z ð x Þ i i i¼1 n Xn 2 ðzðxi Þ−Z * ðxi Þ i¼1 NSE ¼ 1− X n ½ðzðxi Þ−O2 i¼1 Xn ðzðxi Þ−Z * ðxi Þ i¼1 PBIAS ¼ 100: Xn ðzðxi Þ i¼1
ð10Þ ð11Þ ð12Þ
ð13Þ
where z(xi), z*(xi), and O are the measured, interpolated, and mean of the observed values of water quality index of the ith well, respectively, while n is the sample size (Xie et al. 2011). Nash–Sutcliffe efficiencies can range from −∞ to 1. An efficiency of 1 (E=1) corresponds to a perfect match of modeled discharge to the observed data. An efficiency of 0 (E=0) indicates that the model predictions are as accurate as the mean of the observed data, whereas an efficiency less than zero (E< 0) occurs when the observed mean is a better predictor than the model or, in other words, when the residual variance (described by the numerator in the expression above) is larger than the data variance (described by the denominator). Essentially, the closer the model efficiency is to 1, the more accurate the model is (Wagner et al. 2012). Standard deviation (SD) and coefficient of variation (CV) can be used to evaluate uncertainties of different interpolation methods in a certain point. SD is one of the most commonly used indicators of dispersion tendency, while CV is a normalized indicator of dispersion tendency. Smaller SD and CV values indicate less uncertainty in a certain point using different methods. 1Xn x i¼1 i n vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 uX n u x −X t i i¼1 SD ¼ n−1 vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 uX n u x −X t i i¼1 CV ¼ =X n−1
X ¼
ð14Þ
ð15Þ
ð16Þ
where xi is the value in a random point of the ith interpolation method, e.g., n=6. Monte Carlo simulation method was used to analyze the difference among six interpolation methods by setting certain random points in the study area. Different interpolation results were extracted by these random points and analyzed in ArcGIS. SD and CV were used to evaluate the uncertainties of the six different interpolation methods in each random point (Liu et al. 2014).
Assessment of the effect of interpolation method on groundwater contamination Cross-validation only validates the prediction accuracy at a sample site and cannot reflect the spatial difference of interpolation methods. Therefore, in this study, we used the raster analysis function of ESRI ArcGIS to compare the area and spatial differences of contaminated areas estimated by different interpolation methods. In order to analyze the effect of model parameters on contamination assessment, the weighting power of IDW using four levels including 1–4 was applied. For OK, it was required that the distribution of samples to be normal; otherwise, a suitable transformation should be applied to the sample data. Kolmogorov–Smirnov test was used to test the normality of the distributions of water quality index (pIDW2 >OK>IDW3 >IDW4, while the rank of interpolation methods for which the contaminated area was underestimated (level 4) was OK>IDW1> EBK > IDW2 > IDW3 > IDW4, whereas the rank of the methods for which the local minima (level 1) was overestimated was EBK > IDW1 > OK > IDW2 > IDW3 > IDW4. Therefore, it can be concluded that EBK is a better method in comparison with OK because of its less crossvalidation error and its better smoothing effect. Many studies have shown that model parameters are the most important factors for obtaining accurate simulations (Liu et al. 2014; Zhang et al. 2014). In general, IDW needs less input parameters and it is easy to use as well. In contrast, OK needs more parameters and it also has more complex
2767
calculating process (Liu et al. 2014). As the semivariance function fitting is subjective, so different researchers may reach different conclusions (Xie et al. 2011). There is no doubt that OK has the strongest ability to predict the overall trend of groundwater contamination (Xie et al. 2011). It should be noted that the ArcGIS implementation of kriging family of methods, which was adopted in this manuscript, is not statistically rigorous as that of other implementation techniques (e.g., in R packages), but nowadays, many researchers use ArcGIS software for interpolation, so we compared EBK with OK in ArcGIS package. It is recommended that other researchers should be careful on the application of our conclusion to general kriging methods beyond ArcGIS. However, the redeeming features of this methods made it an appealing method as the parameters of a semivariogram structure in OK are being automatically selected using computer code with an appropriate algorithm to minimize some objective functions, such as RMSE; although, it should be noted that OK is an art and depends on the analyst’s skills and the type of data (Gong et al. 2014). The advantage of EBK is that the manual parameter adjusting is eliminated and it automates the most difficult aspects of building a valid kriging model; therefore, conversion of a manual to an automated process might be one of the reasons for the lower uncertainty of EBK.
Conclusion The comparison of interpolation methods to estimate groundwater contamination showed that the three methods have very similar efficiency and did not generate dramatically different results concerning the WQI data, although EBK was more accurate than other methods, with IDW4 having the biggest estimated error. However, we would recommend the EBK method, especially for smaller datasets that do not have enough samples to fit a semivariogram model. Also, the OK method can be subjective, especially when the semivariogram is fitted by eye, and OK may also be very time-consuming compared with EBK. The comparison between interpolation methods suggests that the importance of the interpolation method depends strongly upon the nature of the local wells’ network. In regions of the study area in which wells are sparse, all interpolation methods converge to a similar, narrow range of predictions. When the nearest well for a given location is relatively far away, different interpolation methods will use the same wells for estimation, and therefore, the results will be similar. In areas of the region where wells were dense, different interpolation methods generated completely different predictions of the spatial distribution of groundwater contamination because the used wells weighted differently by different methods. Besides well density, local variation of WQI values of wells is another main factor affecting the uncertainty of different interpolation. The greater the local variation of
2768
WQI data, the greater the uncertainty of interpolation and vice versa. Interpolation accuracy is a relative concept; the criteria are varying with the purpose of interpolation. In the purpose of identifying the groundwater contaminated areas, in this study, we concluded that the local minima and maxima are likely to be smoothed out by three methods. More importantly, we found this smoothing effect, especially for EBK and OK, occurred more in overestimating of local minima than underestimating of local maxima which could be due to the fact that statistical distribution of data is one of most important factors in spatial interpolation. Although, it should be clear that all the interpolation results have errors. Identification of a region as contaminated using an interpolation method is not as easy as it looks and different factors such as scope of the study area, spatial sampling pattern, and statistical distribution of data and even expert skill may affect the results. Acknowledgments The authors are grateful for the help of Andimeshk Health Network and Iran Ministry of Energy for providing us with water quality data.
References Arslan H (2012) Spatial and temporal mapping of groundwater salinity using ordinary kriging and indicator kriging: the case of Bafra Plain, Turkey. Agric Water Manag 113:57–63. doi:10.1016/j.agwat.2012. 06.015 Carl G, Kühn I (2007) Analyzing spatial autocorrelation in species distributions using Gaussian and logit models. Ecol Model 207:159– 170. doi:10.1016/j.ecolmodel.2007.04.024 Cheng K-S, Lin Y-C, Liou J-J (2008) Rain-gauge network evaluation and augmentation using geostatistics. Hydrol Process 22:2554–2564. doi:10.1002/hyp.6851 Dai F, Zhou Q, Lv Z, Wang X, Liu G (2014) Spatial prediction of soil organic matter content integrating artificial neural network and ordinary kriging in Tibetan Plateau. Ecol Indic 45:184–194. doi:10. 1016/j.ecolind.2014.04.003 Falivene O, Cabrera L, Sáez A (2007) Optimum and robust 3D facies interpolation strategies in a heterogeneous coal zone (Tertiary As Pontes Basin, NW Spain). Int J Coal Geol 71:185–208. doi:10. 1016/j.coal.2006.08.008 Falivene O, Cabrera L, Tolosana-Delgado R, Sáez A (2010) Interpolation algorithm ranking using cross-validation and the role of smoothing effect. A coal zone example. Comput Geosci 36:512–519. doi:10. 1016/j.cageo.2009.09.015 Gong G, Mattevada S, O’Bryant SE (2014) Comparison of the accuracy of kriging and IDW interpolations in estimating groundwater arsenic concentrations in Texas. Environ Res 130:59–69. doi:10.1016/j. envres.2013.12.005 Goovaerts P (2000) Estimation or simulation of soil properties? An optimization problem with conflicting criteria. Geoderma 97:165–186. doi:10.1016/S0016-7061(00)00037-9 Gumiere SJ, Lafond JA, Hallema DW, Périard Y, Caron J, Gallichand J (2014) Mapping soil hydraulic conductivity and matric potential for water management of cranberry: characterisation and spatial interpolation methods. Biosyst Eng 128:29–40. doi:10.1016/j. biosystemseng.2014.09.002 Hattermann FF, Wattenbach M, Krysanova V, Wechsung F (2005) Runoff simulations on the macroscale with the ecohydrological model
Environ Sci Pollut Res (2016) 23:2758–2769 SWIM in the Elbe catchment—validation and uncertainty analysis. Hydrol Process 19:693–714. doi:10.1002/hyp.5625 Heistermann M, Kneis D (2011) Benchmarking quantitative precipitation estimation by conceptual rainfall-runoff modeling. Water Resour Res W06514. doi:10.1029/2010WR009153 Horton RK (1965) An index number system for rating water quality. J Water Contam Control Fed 37:300–306 Hosseini E, Gallichand J, Marcotte D (1994) Theoretical andexperimental performance of spatial interpolation methods for soil salinity analysis. Trans ASAE 37:1799–1807 Jager HI, King AW, Schumaker NH, Ashwood TL, Jackson BL (2005) Spatial uncertainty analysis of population models. Ecol Model 185: 13–27. doi:10.1016/j.ecolmodel.2004.10.016 Joseph J, Sharif HO, Sunil T, Alamgir H (2013) Application of validation data for assessing spatial interpolation methods for 8-h ozone or other sparsely monitored constituents. Environ Pollut 178:411– 418. doi:10.1016/j.envpol.2013.03.035 Journel A, Kyriakidis P, Mao S (2000) Correcting the smoothing effect of estimators: a spectral postprocessor. Math Geol 32:787–813. doi:10. 1023/A:1007544406740 Keblouti M, Ouerdachi L, Boutaghane H (2012) Spatial interpolation of annual precipitation in Annaba-Algeria—comparison and evaluation of methods. Energy Procedia 18:468–475. doi:10.1016/j. egypro.2012.05.058 Kravchenko AN (2003) Influence of spatial structure on accuracy of interpolation methods. Soil Sci Soc Am J 67:1564–1571. doi:10. 2136/sssaj2003.1564 Krivoruchko K (2012a) Empirical Bayesian Kriging. ArcUser Fall 2012, pp. 6–10. Also available online at http://www.esri.com/news/ arcuser/1012/empirical-byesian-kriging.html Krivoruchko K (2012b) Modeling Contamination Using Empirical Bayesian Kriging. ArcUser Fall 2012 http://www.esri.com/news/ arcuser/1012/modeling-contamination-using-empirical-bayesiankriging.html Kumar V (2007) Optimal contour mapping of groundwater levels using universal kriging—a case study. Hydrol Sci J 52:1038–1050 Lee J-J, Jang C-S, Wang S-W, Liu C-W (2007) Evaluation of potential health risk of arsenic-affected groundwater using indicator kriging and dose response model. Sci Total Environ 384:151–162. doi:10. 1016/j.scitotenv.2007.06.021 Liu C-W, Jang C-S, Liao C-M (2004) Evaluation of arsenic contamination potential using indicator kriging in the Yun-Lin aquifer (Taiwan). Sci Total Environ 321:173–188. doi:10.1016/j.scitotenv. 2003.09.002 Liu R, Chen Y, Sun C, Zhang P, Wang J, Yu W, Shen Z (2014) Uncertainty analysis of total phosphorus spatial-temporal variations in the Yangtze River Estuary using different interpolation methods. Mar Pollut Bull 86:68–75. doi:10.1016/j.marpolbul.2014.07.041 Lloyd CD (2005) Assessing the effect of integrating elevation data into the estimation of monthly precipitation in Great Britain. J Hydrol 308:128–150. doi:10.1016/j.jhydrol.2004.10.026 Lu GY, Wong DW (2008) An adaptive inverse-distance weighting spatial interpolation technique. Comput Geosci 34:1044–1055. doi:10. 1016/j.cageo.2007.07.010 Mirzaei R, Ghorbani H, Hafezi Moghaddas N, Martín JAR (2014) Ecological risk of heavy metal hotspots in topsoils in the Province of Golestan, Iran. J Geochem Explor 147:268–276. doi:10.1016/j. gexplo.2014.06.011 Nabi Bidhendi GR, Karbassi AR, Nasrabadi T, Hoveidi H (2007) Influence of copper mine on surface water quality. Int J Environ Sci Tech 4:85–91. doi:10.1007/BF03325965 Nampak H, Pradhan B, Manap MA (2014) Application of GIS based data driven evidential belief function model to predict groundwater potential zonation. J Hydrol 513:283–300. doi:10.1016/j.jhydrol.2014. 02.053
Environ Sci Pollut Res (2016) 23:2758–2769 Nas B, Berktay A (2010) Groundwater quality mapping in urban groundwater using GIS. Environ Monit Assess 160:215–227. doi:10.1007/ s10661-008-0689-4 Neshat A, Pradhan B, Javadi S (2015) Risk assessment of groundwater pollution using Monte Carlo approach in an agricultural region: an example from Kerman Plain, Iran. Comput Environ Urban Syst 50: 66–73. doi:10.1016/j.compenvurbsys.2014.11.004 Pan H, Huang W-Q (2003) Influence of uncertainty in delimitation of seismic statistical zone on results of PSHA. Acta Seimol Sin 16: 213–218. doi:10.1007/s11589-003-0024-7 Patel S, Desai KK (2006) Studies on water quality index of some villages of Suratdistrict. Poll Res 25(2):281–285 Plouffe CCF, Robertson C, Chandrapala L (2015) Comparing interpolation techniques for monthly rainfall mapping using multiple evaluation criteria and auxiliary data sources: a case study of Sri Lanka. Environ Model Softw 65:57–71. doi:10.1016/j.envsoft.2015.01.011 Rupal M, Tannushree B, Sukalyan C (2012) Quality characterization of groundwater using water quality index in Surat City, Gujarat, India. Int Res J Environ Sci 1(4):14–23.
2769 Wagner PD, Fiener P, Wilken F, Kumar S, Schneider K (2012) Comparison and evaluation of spatial interpolation schemes for daily rainfall in data scarce regions. J Hydrol 464–465:388–400. doi: 10.1016/j.jhydrol.2012.07.026 Wong DW, Yuan L, Perlin SA (2004) Comparison of spatial interpolation methods for the estimation of air quality data. J Expo Anal Environ Epidemiol 14:404–415. doi:10.1038/sj.jea.7500338 Xie Y, Chen T-B, Lei M, Yang J, Guo Q-J, Song B, Zhou X-Y (2011) Spatial distribution of soil heavy metal pollution estimated by different interpolation methods: accuracy and uncertainty analysis. Chemosphere 82:468–476. doi:10.1016/j. chemosphere.2010.09.053 Yasrebi J, Saffari M, Fathi H, Karimian N, Moazallahi M, Gazni R (2009) Evaluation and comparison of ordinary kriging and inverse distance weighting methods for prediction of spatial variability of some soil chemical parameters. Res J Biol Sci 4:93–102 Zhang P, Liu R, Bao Y, Wang J, Yu W, Shen Z (2014) Uncertainty of SWAT model at different DEM resolutions in a large mountainous watershed. Water Res 53:132–144. doi:10.1016/j.watres. 2014.01.018
RESEARCH ARTICLE
Comparison of interpolation methods for the estimation of groundwater contamination in Andimeshk-Shush Plain, Southwest of Iran Rouhollah Mirzaei 1 & Mohamad Sakizadeh 2
Received: 7 August 2015 / Accepted: 27 September 2015 / Published online: 7 October 2015 # Springer-Verlag Berlin Heidelberg 2015
Abstract Selection of appropriate interpolation methods for the conversion of discrete samples into continuous maps is a controversial issue in the environmental researches. The main objective of this study was to analyze the suitability of three interpolation methods for the discrimination of groundwater with respect to the water quality index (WQI). The groundwater quality data consisted of 17 variables associated with 65 wells located in Andimeshk-Shush Plain. Three spatial interpolation methods including ordinary kriging (OK), empirical Bayesian kriging (EBK), and inverse distance weighting (IDW) were utilized for modeling the groundwater contamination. In addition, different cross-validation indicators were applied to assess the performance of different interpolation methods. The results showed that the performance differed slightly among different methods, although the best performed interpolation method in this study was the empirical Bayesian kriging. Among the interpolation methods, IDW with weighting power of 4 estimated the most contaminated area, while OK estimated the lowest contaminated area. The weighting power of IDW had a significant influence on the estimation, meaning that the estimated contaminated area was increased when a greater weighting power was selected. The subtraction results indicated that there are slightly spatial differences among the contamination assessment results. Results of both standard deviation (SD) and coefficient of variation Responsible editor: Marcus Schulz * Rouhollah Mirzaei [email protected] 1
Department of Environment, Faculty of Natural Resources and Earth Sciences, University of Kashan, Kashan, Iran
2
Department of Environmental Sciences, Faculty of Sciences, Shahid Rajaee Teacher Training University, Tehran, Iran
(CV) also showed that uncertainty was highest in the southern part of the study area, where the distribution of wells were more intensive than that of the northern part. Keywords Groundwater contamination . Uncertainty . Inverse distance weighting . Empirical Bayesian kriging . Iran
Introduction Groundwater is a precious water resource which is used in different fields such as drinking, irrigation, and industrial processes (Nampak et al. 2014; Neshat et al. 2015). As a convenient substitute, groundwater also provides additional water supply to compensate for insufficient surface water resources. As such, groundwater should be protected from contamination by implementing effective water resource planning and management (Neshat et al. 2015). Iran is one of 27 countries that are likely to face increasing water shortage crisis by 2025 unless action is taken to reduce current water consumption (Nabi Bidhendi et al. 2007). Most of the groundwater resources in Iran are used as drinking water besides agricultural and industrial purposes. So, the quality of groundwater is of great importance. Low quality water supply can cause health problems; therefore, determining the quality of water is an important issue. The estimation of spatial distribution of contaminated groundwater is very important in the health risk assessment (Lee et al. 2007). However, only a small proportion of in situ data can be analyzed in a field investigation, due to time and cost constraints. Consequently, sparse measured data contain considerable uncertainty (Liu et al. 2004). Therefore, mapping the spatial distribution of groundwater contamination requires spatial interpolation methods. Surface interpolation methods in a geographic information system (GIS) are very powerful
Environ Sci Pollut Res (2016) 23:2758–2769
tools for predicting surface values (Gong et al. 2014). The discipline of geostatistics provides very useful techniques for handling spatially distributed data such as groundwater contamination (Arslan 2012; Nas and Berktay 2010). By identifying spatial patterns and interpolating values at unsampled locations, interpolation analysis can play a vital role in the sustainable management of groundwater systems by providing estimated input parameters at regular grid points from measurements taken at random locations (Arslan 2012; Kumar 2007). Many authors have emphasized the role of geostatistics in the management and sustainability of regional water resources (Arslan 2012). The most commonly used interpolation methods are inverse distance weighted, kriging, and cokriging interpolations (Gong et al. 2014). However, it has been inconsistent concerning a potential superiority of one interpolation method over another (Gong et al. 2014). Consequently, interpolation methods such as inverse distance weighting (IDW) and kriging have been extensively used in groundwater investigations and contamination mapping (Arslan 2012; Gong et al. 2014). Interpolation accuracy is related to the precise definition of the contaminated area and its boundaries. Consequently, this directly affects the accuracy of contamination assessment. There are a lot of studies on the performance of the spatial interpolation methods for different objectives, but the results are not conclusive. Through interpolation of different data set, some of the authors found that the kriging method has the lowest prediction bias among other methods (Joseph et al. 2013; Liu et al. 2014; Yasrebi et al. 2009), while some other authors showed that IDW was the best method (Gong et al. 2014; Keblouti et al. 2012). Moreover, there are other studies in which different interpolation methods had outperformed such as TPS (Gumiere et al. 2014), ANN kriging (Dai et al. 2014), and Bayesian kriging (Plouffe et al. 2015). Interpolation methods all have a smoothing effect, which underestimates the local high values and overestimates the local low values (Journel et al. 2000; Xie et al. 2011). This smoothing effect leads to bias (e.g., underestimation or overestimation) in contamination assessment and has an effect on relevant environmental decision-making (Goovaerts 2000; Xie et al. 2011). It is essential to minimize this bias in a contamination assessment caused by interpolation methods and to understand the uncertainty of groundwater contamination assessment introduced by interpolation error and differences in contamination assessment between various interpolation methods (Xie et al. 2011). Considering the importance of estimating groundwater contamination, this study was undertaken to compare the performance of three interpolation methods, including IDW, ordinary kriging (OK), and emperical Bayesian kriging (EBK), in the estimation of groundwater quality in Andimeshk-Shush Plain using the leaveone-out cross-validation.
2759
Material and methods Study area Andimeshk-Shush Plain covers an area of approximately 1100 km2 in Khuzestan, Iran (Fig. 1). This plain is located in the north of Khuzestan (32° 00′–32° 35′ N; 48° 10′–48° 25′ E). The local economy depends largely upon farming. Tourism and manufacturing also contribute to the area’s economy. Andimeshk-Shush Plain is comprised of a succession of Dez and Balarood seasonal river deposits interspersed with minor silt in Andimeshk Plain in addition to sand and gravel deposits of Dez and Kharkhe Rivers interspersed with major clayey silt lenses in Shush Plain. Various lithological units range from Cenozoic (Pliocene) to Quaternary Periodare are located within the Andimeshk-Shush Plain. Quaternary-age deposits consist of alluvium, which have been made up of loose, interlayer clay, silt, sand, and gravel. The thickness of the alluvium is about 200–300 m. Andimeshk-Shush aquifer is the primary source of groundwater, supplying approximately 100 % of the total drinking water for about 180,000 people residing in this area. Annual precipitation is approximately 270 mm from which more than 80 % occurs during the December–April period. In comparison, annual potential evaporation is about 1670 mm that is six times higher than that of the annual precipitation. Farms occupy over 70 % of the study area, in which the main agricultural crops are wheat, corn, and sugarcane. Nearly more than 75 % of farmlands in the study area are irrigated by the surface water from Dez irrigation network and the rest (less than 25 %) are irrigated with groundwater. Groundwater quality parameters and water quality index calculation The groundwater quality data used in this study included 17 variables (EC, TDS, turbidity, pH, total hardness, Ca, Mg, sulfate, phosphate, nitrate, nitrite, fluoride, chloride, Fe, Mn, Cu, Cr(VI)) which have been sampled periodically by Andimeshk Health Network and Iran Ministry of Energy from 65 wells for 8 years between 2006 and 2013. Minimum, maximum, average, and standard deviation of these parameters have been given in Table 1. To calculate the water quality index (WQI), the method proposed by Horton (1965) and followed by many researchers Patel and Desai (2006) was used. For computing WQI, three steps were followed as explained by Rupal et al. (2012). We used annual mean of each parameter for WQI calculations. In the first step, a weight (wi) was assigned to each parameter according to its relative importance in the overall quality of water for drinking purposes. The highest weights were assigned to parameters with the most detrimental effect on the health of groundwater consumers, while the parameters with the least issue in
2760
Environ Sci Pollut Res (2016) 23:2758–2769
Fig. 1 A sketch map of the study includes the distribution of sampling wells in Andimeshk-Shush Plain
this regards took the lowest weights. In the second step, the relative weights (Wi) were computed according to this equation: Xn w ð1Þ W i ¼ wi = i¼1 i where Wi is the relative weight, wi is the weight of each parameter, and n is the number of parameters. In the third step, a quality rating scale (qi) for each parameter was assigned by dividing its concentration in each water sample by its respective standard and the result multiplied by 100: qi ¼ C i =S i 100
ð2Þ
where Ci is the concentration of each chemical parameter in each water sample, and Si is the associated drinking water standard. The final WQI was determined by the product of Wi and qi as the following: WQI ¼ W i qi
ð3Þ
Inverse distance weighting Inverse distance weighting is based on the premise that the predictions are a linear combination of available data. In this method, the interpolating function is as follows: Xn w i zi Z ðxÞ ¼ Xi¼1 ð4Þ n w i i¼1 In which wi¼di−u where Z(x) is the predicted value at an interpolated point, whereas Zi is the amount at a known point. n is the total number of known points used in interpolation, di is the distance between point i and the prediction point, wi is the weight assigned to point i. Higher weighting values are assigned to those points which are closer to the interpolated point. As the distance increases, the weight decreases and u is the weighting power that impose the amount of weight decrease with respect to the increase in distance (Xie et al. 2011).
Interpolation methods Ordinary kriging In this study, the most widely used interpolation methods, inverse distance weighting, ordinary kriging, and empirical Bayesian kriging were evaluated.
Kriging is a linear estimator meaning that the estimate of the unknown value is a linear combination of the known data
Environ Sci Pollut Res (2016) 23:2758–2769
2761
Table 1 Descriptive statistics of groundwater quality variables in Andimeshk-Shush Plain along with the associated WQIs Groundwater quality variables
Mean (mg/l)
Min (mg/l)
Max (mg/l)
SD
pH
7.91 0.46 37.69 18.70 240.16 77.06 596.88
6 0.02 0 0 0 8 8
8.13 0.95 102.32 85 1125 213 36,545
2.23 0.189 0.189 17.48 203.56 34.37 4528.60
858.69 435.74
121 60
2367 1185
469.21 229.91
F− Cl− NO3− SO42− Ca2+ Mg2+ EC TDS Turbidity
1.77
0
21
3.40
Total hardness Phosphate
340.68 0.16
73 0.034
1200 0.47
172.53 0.073
Nitrite Fe Mn
0.07 0.09 0.2
0.003 0.01 0
0.46 0.27 2.66
0.107 0.077 0.410
Cu Cr(VI) WQI
0.42 0.5 76.61
0.015 0.004 55
2.90 0.076 90.96
0.646 0.015 8.52
values (Xie et al. 2011). The aim of kriging is to estimate the value of a random function, z, at one or more unsampled points or over larger blocks, from more or less sparse sample data on a given support, say z(x1), z(x2),…z(xN), at x1, x2,…xN. This can be shown by Xn ð5Þ w Z xj z * ð x0 Þ ¼ i¼1 i where wj are the weights assigned to the known value of z(xj) and z*(x0) is the estimated value. To ensure that the estimate is unbiased, weights are made to sum to 1 (Xie et al. 2011): Xn w ¼1 ð6Þ i¼1 i
Emperical Bayesian kriging Empirical Bayesian kriging (EBK) is a geostatistical interpolation method that automates the most difficult aspects of building a valid kriging model. In addition to accounting for the uncertainty in the underlying semivariogram parameters, the other main redeeming feature of EBK is that despite the common applied Geostatistical Analyst in ArcGIS10.2, the parameters in the new developed EBK are automatically optimized through a subsetting and simulation process which is implemented by estimating a lot of semivariogram models
instead of a single semivariogram. The prediction in unknown locations in common kriging methods is done through calculation of semivariogram with respect to the known data locations resulting in the underestimation of the standard error of the prediction due to overlooking the uncertainty of semivariogram. On the contrary, EBK uses an intrinsic random function as the kriging model despite the other kriging methods. The other main difference of EBK with that of the other kriging models is that EBK does not assume a tendency toward an overall mean; thus, there is the same chance for large deviations to get larger or get smaller. The following process is followed in EBK. (1) Using the available data, a semivariogram model is estimated. (2) Given this semivariogram, a new value is simulated at each of the input data location. (3) With respect to the simulated data, a new semivariogram model is estimated accordingly. The calculation of a weight for the latest semivariogram according to Bayes’ rule is the next step in this field. The semivariogram estimated in step 1 is used to simulate a new set of values at the input location d u r i n g t h e r e pe t i t i o n o f s t e p s 2 an d 3 . A ne w semivariogram model and its weight are produced given the simulated data. During this process, the predictions and their respective standard errors are produced at the unsampled locations. This process finally creates a spectrum of semivariograms (Krivoruchko 2012a, b). There are two base distributions available in EBK with respect to the utilized multiplicative skewing normal score in this method: empirical and log empirical. Since log empirical only accept positive data values and their predictions are also positive, so it is a good option for water quality indices that are just positive scores. Due to the application of log transformation on our water quality indices, exponential model was applied to interpolate WQI in this study. Comparison of the accuracy of the interpolation methods Cross-validation and validation with an independent data set are the commonly used methods for comparing the interpolation methods. Because the sample size was limited, a leave-one-out cross-validation (LOOCV) was used to estimate water quality index of each well. With LOOCV, one data point (e.g., water quality index in one well) is removed from the dataset, and its value is estimated by the remaining known values. This process is repeated until water quality index values from all the wells are estimated (Xie et al. 2011). In order to assess the performance of different interpolation methods the mean error (ME), mean absolute error (MAE), mean relative errors (MRE), mean squared error (MSE), root mean squared errors (RMSE), Nash–Sutcliffe efficiency
2762
Environ Sci Pollut Res (2016) 23:2758–2769
(NSE), and percentage bias (PBIAS) of the interpolated values were calculated with the following equations: Xn Z * ðxi Þ−Z ðxi Þ i¼1 ð7Þ ME ¼ n Xn Z * ðxi Þ−Z ðxi Þ i¼1 MAE ¼ ð8Þ n 1 X n Z * ðxi Þ−Z ðxi Þ ð9Þ MRE ¼ i¼1 n Z ðx Þ X
MSE ¼
i
Z ðxi Þ−zðxi Þ *
2
n rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X n 2 1 * RMSE ¼ Z ð x Þ−Z ð x Þ i i i¼1 n Xn 2 ðzðxi Þ−Z * ðxi Þ i¼1 NSE ¼ 1− X n ½ðzðxi Þ−O2 i¼1 Xn ðzðxi Þ−Z * ðxi Þ i¼1 PBIAS ¼ 100: Xn ðzðxi Þ i¼1
ð10Þ ð11Þ ð12Þ
ð13Þ
where z(xi), z*(xi), and O are the measured, interpolated, and mean of the observed values of water quality index of the ith well, respectively, while n is the sample size (Xie et al. 2011). Nash–Sutcliffe efficiencies can range from −∞ to 1. An efficiency of 1 (E=1) corresponds to a perfect match of modeled discharge to the observed data. An efficiency of 0 (E=0) indicates that the model predictions are as accurate as the mean of the observed data, whereas an efficiency less than zero (E< 0) occurs when the observed mean is a better predictor than the model or, in other words, when the residual variance (described by the numerator in the expression above) is larger than the data variance (described by the denominator). Essentially, the closer the model efficiency is to 1, the more accurate the model is (Wagner et al. 2012). Standard deviation (SD) and coefficient of variation (CV) can be used to evaluate uncertainties of different interpolation methods in a certain point. SD is one of the most commonly used indicators of dispersion tendency, while CV is a normalized indicator of dispersion tendency. Smaller SD and CV values indicate less uncertainty in a certain point using different methods. 1Xn x i¼1 i n vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 uX n u x −X t i i¼1 SD ¼ n−1 vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 uX n u x −X t i i¼1 CV ¼ =X n−1
X ¼
ð14Þ
ð15Þ
ð16Þ
where xi is the value in a random point of the ith interpolation method, e.g., n=6. Monte Carlo simulation method was used to analyze the difference among six interpolation methods by setting certain random points in the study area. Different interpolation results were extracted by these random points and analyzed in ArcGIS. SD and CV were used to evaluate the uncertainties of the six different interpolation methods in each random point (Liu et al. 2014).
Assessment of the effect of interpolation method on groundwater contamination Cross-validation only validates the prediction accuracy at a sample site and cannot reflect the spatial difference of interpolation methods. Therefore, in this study, we used the raster analysis function of ESRI ArcGIS to compare the area and spatial differences of contaminated areas estimated by different interpolation methods. In order to analyze the effect of model parameters on contamination assessment, the weighting power of IDW using four levels including 1–4 was applied. For OK, it was required that the distribution of samples to be normal; otherwise, a suitable transformation should be applied to the sample data. Kolmogorov–Smirnov test was used to test the normality of the distributions of water quality index (pIDW2 >OK>IDW3 >IDW4, while the rank of interpolation methods for which the contaminated area was underestimated (level 4) was OK>IDW1> EBK > IDW2 > IDW3 > IDW4, whereas the rank of the methods for which the local minima (level 1) was overestimated was EBK > IDW1 > OK > IDW2 > IDW3 > IDW4. Therefore, it can be concluded that EBK is a better method in comparison with OK because of its less crossvalidation error and its better smoothing effect. Many studies have shown that model parameters are the most important factors for obtaining accurate simulations (Liu et al. 2014; Zhang et al. 2014). In general, IDW needs less input parameters and it is easy to use as well. In contrast, OK needs more parameters and it also has more complex
2767
calculating process (Liu et al. 2014). As the semivariance function fitting is subjective, so different researchers may reach different conclusions (Xie et al. 2011). There is no doubt that OK has the strongest ability to predict the overall trend of groundwater contamination (Xie et al. 2011). It should be noted that the ArcGIS implementation of kriging family of methods, which was adopted in this manuscript, is not statistically rigorous as that of other implementation techniques (e.g., in R packages), but nowadays, many researchers use ArcGIS software for interpolation, so we compared EBK with OK in ArcGIS package. It is recommended that other researchers should be careful on the application of our conclusion to general kriging methods beyond ArcGIS. However, the redeeming features of this methods made it an appealing method as the parameters of a semivariogram structure in OK are being automatically selected using computer code with an appropriate algorithm to minimize some objective functions, such as RMSE; although, it should be noted that OK is an art and depends on the analyst’s skills and the type of data (Gong et al. 2014). The advantage of EBK is that the manual parameter adjusting is eliminated and it automates the most difficult aspects of building a valid kriging model; therefore, conversion of a manual to an automated process might be one of the reasons for the lower uncertainty of EBK.
Conclusion The comparison of interpolation methods to estimate groundwater contamination showed that the three methods have very similar efficiency and did not generate dramatically different results concerning the WQI data, although EBK was more accurate than other methods, with IDW4 having the biggest estimated error. However, we would recommend the EBK method, especially for smaller datasets that do not have enough samples to fit a semivariogram model. Also, the OK method can be subjective, especially when the semivariogram is fitted by eye, and OK may also be very time-consuming compared with EBK. The comparison between interpolation methods suggests that the importance of the interpolation method depends strongly upon the nature of the local wells’ network. In regions of the study area in which wells are sparse, all interpolation methods converge to a similar, narrow range of predictions. When the nearest well for a given location is relatively far away, different interpolation methods will use the same wells for estimation, and therefore, the results will be similar. In areas of the region where wells were dense, different interpolation methods generated completely different predictions of the spatial distribution of groundwater contamination because the used wells weighted differently by different methods. Besides well density, local variation of WQI values of wells is another main factor affecting the uncertainty of different interpolation. The greater the local variation of
2768
WQI data, the greater the uncertainty of interpolation and vice versa. Interpolation accuracy is a relative concept; the criteria are varying with the purpose of interpolation. In the purpose of identifying the groundwater contaminated areas, in this study, we concluded that the local minima and maxima are likely to be smoothed out by three methods. More importantly, we found this smoothing effect, especially for EBK and OK, occurred more in overestimating of local minima than underestimating of local maxima which could be due to the fact that statistical distribution of data is one of most important factors in spatial interpolation. Although, it should be clear that all the interpolation results have errors. Identification of a region as contaminated using an interpolation method is not as easy as it looks and different factors such as scope of the study area, spatial sampling pattern, and statistical distribution of data and even expert skill may affect the results. Acknowledgments The authors are grateful for the help of Andimeshk Health Network and Iran Ministry of Energy for providing us with water quality data.
References Arslan H (2012) Spatial and temporal mapping of groundwater salinity using ordinary kriging and indicator kriging: the case of Bafra Plain, Turkey. Agric Water Manag 113:57–63. doi:10.1016/j.agwat.2012. 06.015 Carl G, Kühn I (2007) Analyzing spatial autocorrelation in species distributions using Gaussian and logit models. Ecol Model 207:159– 170. doi:10.1016/j.ecolmodel.2007.04.024 Cheng K-S, Lin Y-C, Liou J-J (2008) Rain-gauge network evaluation and augmentation using geostatistics. Hydrol Process 22:2554–2564. doi:10.1002/hyp.6851 Dai F, Zhou Q, Lv Z, Wang X, Liu G (2014) Spatial prediction of soil organic matter content integrating artificial neural network and ordinary kriging in Tibetan Plateau. Ecol Indic 45:184–194. doi:10. 1016/j.ecolind.2014.04.003 Falivene O, Cabrera L, Sáez A (2007) Optimum and robust 3D facies interpolation strategies in a heterogeneous coal zone (Tertiary As Pontes Basin, NW Spain). Int J Coal Geol 71:185–208. doi:10. 1016/j.coal.2006.08.008 Falivene O, Cabrera L, Tolosana-Delgado R, Sáez A (2010) Interpolation algorithm ranking using cross-validation and the role of smoothing effect. A coal zone example. Comput Geosci 36:512–519. doi:10. 1016/j.cageo.2009.09.015 Gong G, Mattevada S, O’Bryant SE (2014) Comparison of the accuracy of kriging and IDW interpolations in estimating groundwater arsenic concentrations in Texas. Environ Res 130:59–69. doi:10.1016/j. envres.2013.12.005 Goovaerts P (2000) Estimation or simulation of soil properties? An optimization problem with conflicting criteria. Geoderma 97:165–186. doi:10.1016/S0016-7061(00)00037-9 Gumiere SJ, Lafond JA, Hallema DW, Périard Y, Caron J, Gallichand J (2014) Mapping soil hydraulic conductivity and matric potential for water management of cranberry: characterisation and spatial interpolation methods. Biosyst Eng 128:29–40. doi:10.1016/j. biosystemseng.2014.09.002 Hattermann FF, Wattenbach M, Krysanova V, Wechsung F (2005) Runoff simulations on the macroscale with the ecohydrological model
Environ Sci Pollut Res (2016) 23:2758–2769 SWIM in the Elbe catchment—validation and uncertainty analysis. Hydrol Process 19:693–714. doi:10.1002/hyp.5625 Heistermann M, Kneis D (2011) Benchmarking quantitative precipitation estimation by conceptual rainfall-runoff modeling. Water Resour Res W06514. doi:10.1029/2010WR009153 Horton RK (1965) An index number system for rating water quality. J Water Contam Control Fed 37:300–306 Hosseini E, Gallichand J, Marcotte D (1994) Theoretical andexperimental performance of spatial interpolation methods for soil salinity analysis. Trans ASAE 37:1799–1807 Jager HI, King AW, Schumaker NH, Ashwood TL, Jackson BL (2005) Spatial uncertainty analysis of population models. Ecol Model 185: 13–27. doi:10.1016/j.ecolmodel.2004.10.016 Joseph J, Sharif HO, Sunil T, Alamgir H (2013) Application of validation data for assessing spatial interpolation methods for 8-h ozone or other sparsely monitored constituents. Environ Pollut 178:411– 418. doi:10.1016/j.envpol.2013.03.035 Journel A, Kyriakidis P, Mao S (2000) Correcting the smoothing effect of estimators: a spectral postprocessor. Math Geol 32:787–813. doi:10. 1023/A:1007544406740 Keblouti M, Ouerdachi L, Boutaghane H (2012) Spatial interpolation of annual precipitation in Annaba-Algeria—comparison and evaluation of methods. Energy Procedia 18:468–475. doi:10.1016/j. egypro.2012.05.058 Kravchenko AN (2003) Influence of spatial structure on accuracy of interpolation methods. Soil Sci Soc Am J 67:1564–1571. doi:10. 2136/sssaj2003.1564 Krivoruchko K (2012a) Empirical Bayesian Kriging. ArcUser Fall 2012, pp. 6–10. Also available online at http://www.esri.com/news/ arcuser/1012/empirical-byesian-kriging.html Krivoruchko K (2012b) Modeling Contamination Using Empirical Bayesian Kriging. ArcUser Fall 2012 http://www.esri.com/news/ arcuser/1012/modeling-contamination-using-empirical-bayesiankriging.html Kumar V (2007) Optimal contour mapping of groundwater levels using universal kriging—a case study. Hydrol Sci J 52:1038–1050 Lee J-J, Jang C-S, Wang S-W, Liu C-W (2007) Evaluation of potential health risk of arsenic-affected groundwater using indicator kriging and dose response model. Sci Total Environ 384:151–162. doi:10. 1016/j.scitotenv.2007.06.021 Liu C-W, Jang C-S, Liao C-M (2004) Evaluation of arsenic contamination potential using indicator kriging in the Yun-Lin aquifer (Taiwan). Sci Total Environ 321:173–188. doi:10.1016/j.scitotenv. 2003.09.002 Liu R, Chen Y, Sun C, Zhang P, Wang J, Yu W, Shen Z (2014) Uncertainty analysis of total phosphorus spatial-temporal variations in the Yangtze River Estuary using different interpolation methods. Mar Pollut Bull 86:68–75. doi:10.1016/j.marpolbul.2014.07.041 Lloyd CD (2005) Assessing the effect of integrating elevation data into the estimation of monthly precipitation in Great Britain. J Hydrol 308:128–150. doi:10.1016/j.jhydrol.2004.10.026 Lu GY, Wong DW (2008) An adaptive inverse-distance weighting spatial interpolation technique. Comput Geosci 34:1044–1055. doi:10. 1016/j.cageo.2007.07.010 Mirzaei R, Ghorbani H, Hafezi Moghaddas N, Martín JAR (2014) Ecological risk of heavy metal hotspots in topsoils in the Province of Golestan, Iran. J Geochem Explor 147:268–276. doi:10.1016/j. gexplo.2014.06.011 Nabi Bidhendi GR, Karbassi AR, Nasrabadi T, Hoveidi H (2007) Influence of copper mine on surface water quality. Int J Environ Sci Tech 4:85–91. doi:10.1007/BF03325965 Nampak H, Pradhan B, Manap MA (2014) Application of GIS based data driven evidential belief function model to predict groundwater potential zonation. J Hydrol 513:283–300. doi:10.1016/j.jhydrol.2014. 02.053
Environ Sci Pollut Res (2016) 23:2758–2769 Nas B, Berktay A (2010) Groundwater quality mapping in urban groundwater using GIS. Environ Monit Assess 160:215–227. doi:10.1007/ s10661-008-0689-4 Neshat A, Pradhan B, Javadi S (2015) Risk assessment of groundwater pollution using Monte Carlo approach in an agricultural region: an example from Kerman Plain, Iran. Comput Environ Urban Syst 50: 66–73. doi:10.1016/j.compenvurbsys.2014.11.004 Pan H, Huang W-Q (2003) Influence of uncertainty in delimitation of seismic statistical zone on results of PSHA. Acta Seimol Sin 16: 213–218. doi:10.1007/s11589-003-0024-7 Patel S, Desai KK (2006) Studies on water quality index of some villages of Suratdistrict. Poll Res 25(2):281–285 Plouffe CCF, Robertson C, Chandrapala L (2015) Comparing interpolation techniques for monthly rainfall mapping using multiple evaluation criteria and auxiliary data sources: a case study of Sri Lanka. Environ Model Softw 65:57–71. doi:10.1016/j.envsoft.2015.01.011 Rupal M, Tannushree B, Sukalyan C (2012) Quality characterization of groundwater using water quality index in Surat City, Gujarat, India. Int Res J Environ Sci 1(4):14–23.
2769 Wagner PD, Fiener P, Wilken F, Kumar S, Schneider K (2012) Comparison and evaluation of spatial interpolation schemes for daily rainfall in data scarce regions. J Hydrol 464–465:388–400. doi: 10.1016/j.jhydrol.2012.07.026 Wong DW, Yuan L, Perlin SA (2004) Comparison of spatial interpolation methods for the estimation of air quality data. J Expo Anal Environ Epidemiol 14:404–415. doi:10.1038/sj.jea.7500338 Xie Y, Chen T-B, Lei M, Yang J, Guo Q-J, Song B, Zhou X-Y (2011) Spatial distribution of soil heavy metal pollution estimated by different interpolation methods: accuracy and uncertainty analysis. Chemosphere 82:468–476. doi:10.1016/j. chemosphere.2010.09.053 Yasrebi J, Saffari M, Fathi H, Karimian N, Moazallahi M, Gazni R (2009) Evaluation and comparison of ordinary kriging and inverse distance weighting methods for prediction of spatial variability of some soil chemical parameters. Res J Biol Sci 4:93–102 Zhang P, Liu R, Bao Y, Wang J, Yu W, Shen Z (2014) Uncertainty of SWAT model at different DEM resolutions in a large mountainous watershed. Water Res 53:132–144. doi:10.1016/j.watres. 2014.01.018
