Study on spatial variation of any groundwater quality parameter is not possible without taking the spatial relationship amongst the parameter values into account. Geostatistical methods are useful tools for conducting such studies. In this study, due to significant spatial and temporal variation of salinity (EC) in Feyz Abad-Mahvelat plain, EC data collected between the years 1386 and 1391 from operational wells in the plain, were analyzed to identify the best Geostatistical method applicable to the EC data in the plain. Although EC values change over time and space, for the spatial analysis an average annual data was considered. Based on the preprocessing of the data, capabilities and properties of different Geostatistical methods, and literature review on the subject, six different Geostatistical methods were selected for the analysis including appropriate types of Ordinary and Universal Kriging and Empirical Bayesian Kriging (EBK). Since there are numerous internal parameters in different methods of Kriging, a systematic multi-step procedure was used to optimize the parameters and to find the best model. To achieve this, 390 combinations of methods and parameters were generated and used in the comparative analysis. First the best models were selected using cross validation and a defined lumped statistical measure that combines two conventional error measures in a single one. Then, sensitivity analysis was performed on the selected best models and also they were evaluated in terms of compliance with the field conditions. The procedure was applied to the first and last year of the available data, separately (year 1386 and year 1391). The results indicated the superiority of the EBK method in terms of both accuracy and performance compared to the other methods. In addition, the generated maps showed the high level of salinity in the northern, southern, and western plain borders and significant salt water intrusion from the western border over time.