A R T I C L E  I N F O		ABSTRACT
ORIGINAL ARTICLE		Introduction: It is generally accepted that groundwater is one of the most vital sources of water for drinking use in cities and rural areas. The water drawn from these sources should be sanitary, have low soluble substances, and be free of any pathogens and microorganisms. Materials and Methods: In this study, 90 wells were sampled with proper dispersion over the study area to achieve suitable estimation accuracy. Results: The assessments made based on 10-year averages of water quality in the studied plain showed that according to the Schuler and Wilcox criteria of water quality for drinking and agricultural use, the northern and southern parts of the plain have unsuitable water quality compared to central parts. Interpolation RMSE value of the Inverse Distance Weighting (IDW) model for SO4-, TDS, TH, Mg2+, Cl-, Ca2+ and HCO3- were 8.46, 2615, 246, 6.8, 38.9, 8.3 1.16 also 8.4, 2628, 750.9, 7.0, 39.8, 8.1, 8.1 (mg L-1) for Kriging, respectively. Conclusion: The cause of low groundwater quality in northern regions is the high rate of SO4-, TH, Cl-, and TDS, which are of the most important determinants of water quality for drinking. The examination of samples in the assessment of water quality for agricultural use clearly showed a higher value of EC compared to SAR.
Article History: Received: 12 August 2019 Accepted: 20 October 2019
*Corresponding Author: Seyed Ali Almodaresi Email: Almodaresi@iauyazd.ac.ir Tel: +989131526455
Keywords: Schuler Diagram, Wilcox Diagram, IDW and Kriging Models, Drinking Water, Ardakan Mehriz plain.

• Evaluation and sampling of 93 wells in the plain and identifying the elements that affect the water quality assessment.
• Kriging interpolation for each of the elements that affect the water quality assessment (obtained from the previous step) according to the relevant contour plot using the inverse distance interpolation method and by an evaluation based on mean squared error.
• Determining the range of water quality by comparing the data obtained from the previous steps using the Schuler diagram, which classifies water into six categories: good, acceptable, moderate, unsuitable, quite unpleasant, and non-potable.
• Zoning the area based on the effective parameters with the help of the Wilcox diagram.

Figure 1: Location of study area of Yazd-Ardakan-Mehriz
 

Data and Method
Geostatistics is one of the best applications of interpolation operation, which has been able to improve the accuracy of estimation through scientific models. Accordingly, geostatistics can be used to convert observations into data. In the GIS environment, geographical data are divided into two types based on the analytical methods: discrete data and continuous data. Discrete data are mostly categorical data, meaning that in essence they have precisely definable boundaries and can be stored as a raster or vector, (e.g. lake, building, or road). Continuous data are continuous in nature, meaning that every point on the ground will have a degree of the property or variable they describe. Liquids belong to the latter group, but also have a direction and can be measured at any point. Because of their continuity, continuous data cannot be measured at all levels and have to be sampled. Interpolation is the estimation of continuous variables in the non-sampled areas within an area with scattered observation points. The output of the interpolation can be used as a plot or layer in the GIS analysis. There are many methods for interpolation of continuous data, but this study used two methods for this purpose: IDW and Kriging models.
Inverse Distance Weighting (IDW)
In IDW model, it is assumed that the effectiveness of the continuous variable decreases with the distance from the unknown point, therefore the distance is used as the weight of the known variable in the estimation of unmeasured points. This method is called Inverse Distance Weighting because as the distance from the unknown point increases, the weight decreases ^{13, 14}. The effect of the intensity of spatial dependence of data can be applied to the formulation by raising the power of the inverse distance. The square of the inverse distance is a popular choice for this purpose. In this method, the target range is transformed into a matrix with cells of equal size. The spatial coordinates of this matrix and the measurement unit are known. For example, a matrix may have 50×50 m² cells. In this grid, the value of the variable is known in some cells and unknown in others. Each unknown value is estimated using the values of neighboring cells within a given radius using the following equation ^{14, 15}.

Where the value is measured at the i^th position and  is the weight measured at the i^th position. S₀ is the position of the estimated value and N is the number of measured or known points. is a function of the distance between points. The dependence of the unknown cells on known cells is adjusted by the power of the inverse distance. The appropriate power ( ) is determined based on the Root Mean Square Prediction Error (RMSPE). The best value for is the one that yields the best estimations of unknown cells, or in the other words minimizes the prediction error. The points on the curve represent different RMSPE values obtained with different  values. The minimum value of RMSPE on the curve gives the best power for the IDW model.
When this power is zero ( ), the impact of the distance is uniform and the unknown value is obtained by averaging the values of neighboring points. But as this power increases, the impact of distance also increases, and closer cells gain a higher weight in the formulation.
In the IDW model, it is typical to use  greater than 1, such as 2, which is why it is also called squared inverse distance. In the IDW model Neighborhood can be defined in two ways. The fixed search radius method assumes a circle around the unknown point and calculates this value based on the points positioned within the circle. In other words, to calculate the unknown point, the distance of each of the points inside the circle will be measured, the measured distances will be raised to second order, and the results will be averaged. The value is actually the weight given to the distances. This means to increase the effect of distance, so the value of is increased. The size of the search radius depends on the distance between the points and the manner of changes in the continuous variable. If the variable changes irregularly, one can use the nearest neighbor’s method instead of the fixed search radius. This method is similar to the previous method, with the difference that interpolation is performed with the minimum number of neighbors. In other words, the neighborhood is defined with the number of neighbors ^{16, 17}.
Kriging Model
The Kriging method is the most important and extensively used method of statistical interpolation. Kriging is an advanced interpolation method most suitable for data with well-defined local trends. The Kriging method operates based on the minimum variance of estimation, and its error is a function of variogram characteristics (spatial structure) ^{18, 19}. As long as variographic analysis and variogram model are sufficiently accurate, the kriging method will be accurate as well. The Kriging model is generally similar to the IDW model ^{20, 21}:

Where  is the value measured at the i^th is position and   is the weight measured at the i^th position. S₀ is the position of the estimated value and N is the number of measured or known points. In the IDW model, is only a function of distance, but in the kriging model, the weight also depends on the spatial structure of the points. Accordingly, the Kriging model is a good interpolation model for geostatistical purposes. The kriging model operates based on the regionalized variable theory, where the value of each point in space is a function of its coordinates. Therefore, the variation of the regionalized variable in space is decomposed into 4 components.

Where m(x) is the function giving the structural element Z at point x,  is the function describing the random component, which is spatially dependent and locally variable, and is the residue or random error, which has a zero mean and constant variance and is spatially independent. In a semi-variogram, every two samples can be paired with each other, as there is a known distance between them. Therefore, for better management of the semivariogram, instead of drawing all pairs, they are grouped based on their distance. For example, all pairs that have a distance of more than 40 meters and less than 50 meters from each other can be categorize in one group. After drawing the semivariograms, a suitable regression model will be fitted into it ^{22, 23}.
Results
Accuracy evaluation method
Different interpolation methods are evaluated using the cross-validation method. In this method, a point will be temporarily removed and will then be estimated based on other points using interpolation. The removed value will then be returned and this process will be repeated for all remaining points. In the end, a table with two columns, one for actual values and the other for estimated values will be created. Having these two values, Mean Absolute Error (MAE) and Mean Bias Error (MBE) of the model can be measured. The closer the MAE and MBE values ​​are to zero, the more accurate is the model. Other measures for assessing the accuracy of interpolation methods are the Root Mean Square Error (RMSE) and the correlation coefficient​​ (R²) between observed and measured values. Greater RMSE values represent a higher error, but in the case of R², greater values indicate higher statistical accuracy ^{24, 25} (Table 1).
These measures are calculated as follows:
               

RMSE of Kriging	RMSE of IDW	Interpolation
8.45	8.46	SO4- (mg L-1)
2628.92	2615.69	Total dissolved solids ( mg L-1)
750.93	746.27	TH (mg L-1)
7.01	6.83	Mg2+ (mg L-1)
39.80	38.90	Cl- (mg L-1)
8.14	8.29	Ca2+ (mg L-1)
8.14	1.16	HCO3- (mg L-1)

Using the samples taken in a 10-year period, the contour plots were obtained for each parameter, and the water quality plots were drawn accordingly. According to Schuler criteria for drinking water quality, it seems that central and southwestern parts of the studied region have good and acceptable water quality, but southeastern and northern parts have a low water quality, which fall into unsuitable and even unpleasant category (Figure 2).

A

B

D

C

E

F

H

G

Water quality assessment for agricultural use
The present study used this classification with the Wilcox criteria. The figures below show the zone plots related to SAR and EC and the results regarding the water quality for irrigation based on the Wilcox diagram and criteria ^26-28 (Figure 3, 4).

B

A

1. Biglari H, Chavoshani A, Javan N, et al. Geochemical study of groundwater conditions with special emphasis on fluoride concentration, Iran. Desalin Water Treat. 2016;57(47):22392-9.
2. Rahmati O, Pourghasemi HR, Melesse AM. Application of GIS-based data driven random forest and maximum entropy models for groundwater potential mapping: a case study at Mehran Region, Iran. Catena. 2016;137:360-72.
3. Izady A, Davary K, Alizadeh A, et al. Groundwater conceptualization and modeling using distributed SWAT-based recharge for the semi-arid agricultural Neishaboor plain, Iran. Hydrogeol. J. 2015;23(1):47-68.
4. Tian Y, Zheng Y, Wu B, et al. Modeling surface water-groundwater interaction in arid and semi-arid regions with intensive agriculture. Environmental Modelling & Software. 2015;63:170-84.
5. Gibbs RJ. Mechanisms controlling world water chemistry. Science. 1970;170(3962):1088-90.
6. Elbeih SF. An overview of integrated remote sensing and GIS for groundwater mapping in Egypt. Ain Shams Eng. J. 2015;6(1):1-15.
7. Van der Lee J, Gehrels J. Modelling of groundwater recharge for a fractured dolomite aquifer under semi-arid conditions.  Recharge of Phreatic Aquifers in (Semi-) Arid Areas: Routledge. 2017: 129-44.
8. Chung S, Venkatramanan S, Kim T, et al. Influence of hydrogeochemical processes and assessment of suitability for groundwater uses in Busan City, Korea. Environment, Development and Sustainability. 2015;17(3):423-41.
9. Nnadi EO, Newman AP, Coupe SJ, et al. Stormwater harvesting for irrigation purposes: An investigation of chemical quality of water recycled in pervious pavement system. J Environ Manage. 2015;147:246-56.
10. Maghami y, Ghazavi R, Vali a. Evaluation of different interpolation methods for water quality zoning using GIS (case study: Abadeh county). J Geog & Environ Planning. 2012;22(42):171-82.
11. Xiao Y, Gu X, Yin S, et al. Geostatistical interpolation model selection based on ArcGIS and spatio-temporal variability analysis of groundwater level in piedmont plains, northwest China. SpringerPlus. 2016;5(1):425.
12. Mirzaei R, Sakizadeh M. Comparison of interpolation methods for the estimation of groundwater contamination in Andimeshk-Shush Plain, Southwest of Iran. Environ Sci Pollut Res. 2016;23(3):2758-69.
13. Varatharajan R, Manogaran G, Priyan M, et al. Visual analysis of geospatial habitat suitability model based on inverse distance weighting with paired comparison analysis. Multimedia tools and applications. 2018; 77(14):17573-93.
14. Yasrebi AB, Afzal P, Wetherelt A, et al. Application of an inverse distance weighted anisotropic method (IDWAM) to estimate elemental distribution in Eastern Kahang Cu-Mo porphyry deposit, Central Iran. International Journal of Mining and Mineral Engineering. 2016;7(4):340-62.
15. Giri S, Qiu Z. Understanding the relationship of land uses and water quality in Twenty First Century: A review. J Environ Manage. 2016;173:41-8.
16. Zhang J, Howard K, Langston C, et al. Multi-Radar Multi-Sensor (MRMS) quantitative precipitation estimation: Initial operating capabilities. Bull Am Meteorol Soc. 2016;97(4): 621-38.
17. Austin PC, Stuart EA. Moving towards best practice when using inverse probability of treatment weighting (IPTW) using the propensity score to estimate causal treatment effects in observational studies. Stat Med. 2015; 34(28): 3661-79.
18. Zhang L, Lu Z, Wang P. Efficient structural reliability analysis method based on advanced Kriging model. Appl Math Modell. 2015;39(2):781-93.
19. Zhao Y, Lu W, Xiao C. A Kriging surrogate model coupled in simulation–optimization approach for identifying release history of groundwater sources. J Contam Hydrol. 2016;185:51-60.
20. Gao Z, Shao X, Jiang P, et al. Parameters optimization of hybrid fiber laser-arc butt welding on 316L stainless steel using Kriging model and GA. Opt Laser Technol. 2016;83:153-62.
21. Mukhopadhyay T, Chakraborty S, Dey S, et al. A critical assessment of Kriging model variants for high-fidelity uncertainty quantification in dynamics of composite shells. Arch Comput Methods Eng. 2017;24(3):495-518.
22. Fan F, Bell K, Hill D, et al, editors. Wind forecasting using kriging and vector auto-regressive models for dynamic line rating studies. Power Tech. 2015.
23. Harman BI, Koseoglu H, Yigit CO. Performance evaluation of IDW, Kriging and multiquadric interpolation methods in producing noise mapping: A case study at the city of Isparta, Turkey. Applied Acoustics. 2016;112: 147-57.
24. Dube T, Mutanga O. Evaluating the utility of the medium-spatial resolution Landsat 8 multispectral sensor in quantifying aboveground biomass in uMgeni catchment, South Africa. J Photogramm Remote Sens. 2015;101:36-46.
25. Funk C, Verdin A, Michaelsen J, et al. A global satellite assisted precipitation climatology. Earth Syst Sci Data Discuss. 2015;8(1):794-9.
26. Alavi N, Zaree E, Hassani M, et al. Water quality assessment and zoning analysis of Dez eastern aquifer by Schuler and Wilcox diagrams and GIS. Desalin Water Treat. 2016;57(50):23686-97.
27. Choramin M, Safaei A, Khajavi S, et al. Analyzing and studding chemical water quality parameters and its changes on the base of Schuler, Wilcox and Piper diagrams (project: Bahamanshir River). WALIA journal. 2015;31:22-7.
28. Fhanduzi F, Pari Zanganeh A, Aamani A, et al. Survey of hydro-geochemical quality and health of groundwater in Ramian, Golestan Province, Iran. Journal of health research in community. 2015;1(3):41-52.