1-Introduction Selection and accuracy of appropriate zoning methods and preparation of maps corresponding to variation of groundwater qualitative characteristics depend upon the region's condition and the available statistics and data. The goal of this research is to determine the most appropriate method of interpolation and spatial analysis of qualitative components such as Cl, EC, SO4, SAR, TDS and TH of groundwater in Mashhad Plain. In this research, at first, the data of 177 wells were selected with respect to their dispersion and accuracy for two consecutive years (1392-1393); then controlling and reconstruction of data were performed. Kolmogorov-Smirnov test showed that the data were not normal and to normalize them their logarithm was obtained. Then using GS+software, the best variogram model was fitted to the spatial structure of data. The results show that for qualitative parameters, the exponential, Semi-variogram and linear models are the best semi variogram models. The precision of groundwater qualitative parameters in the three methods of Kriging, Cokriging and Inverse Distance Weighting (IDW) was assessed and the precision rate results showed that Cokriging method has a higher precision and lower error. The spatial structure of studied characteristics follows the exponential, spHerical and linear model with an influence range of 9656 to 71100001m and the threshold range of 0.965 to 5.65, and the spatial dependency class was in the range of 0-0.87. Finally, the zoning maps of water qualitative parameters were prepared using GIS software. Due to population concentration in the southern parts and overuse of wells and recent droughts, the highest values of qualitative components belong to the south of the study area. According to the geological formations and the type of rocks in the northwest and southeast, the concentration of calcium, potassium, sodium, chlorine, as well as the percentage of evaporation and transpiration and temperature are more evident. 2- Methodology In this research, three methods of interpolation, namely, Weighted distance image, Kriging and Co-kriging are used to predict some quality indicators such as Na, TH, EC, SAR, CL, CA, Mg, and SO4 in 2013-2014. Data have been collected from 177 piezometric wells in the Mashhad Plain. In this study, after data normalization, parameters were interpolated using three geostatistical methods (ordinary kriging, co-kriging, and IDW). In order to evaluate the methods of interpolation, the mutual evaluation test was used. Then, based on the best method of mediation and the use of geographic information system, a map of some parameters and plain zoning maps were prepared. The purpose of this research was to determine the most appropriate interpolation method in order to investigate and spatially analyze the changes in some quality characteristics of underground water in Mashhad plain. After the appropriate model was fitted, according to CO/(CO+C) observations, the strength of the spatial structure in the quality components of the studied water is not very strong. For the zoning of the quality components, the Co-kriging method was introduced as a suitable method, and considering that correlation between different variables is used in this method. It can be said that there is a significant relationship between the components of groundwater (Zahtabian et al., 2019). The results of zoning the quality of the parameters showed two inappropriate ranges in the northwest and east of the region, which about 70% of the exploitation wells are located in these areas. Also, due to the recent droughts, the process of desertification will be accelerated in this region; therefore, more attention should be paid to the proper management of underground water in this region. 4- Discussion & Conclusions According to the obtained results, the best fitted model for CA and MG parameters is linear model; EC parameters are exponential model and NA parameters are Semi-variogram model. Regarding the spatial structure, parameters MG and TH have strong location continuity (less than 0.25) and parameters MG, NA, EC have weak continuity (more than 0.75). According to the obtained results, was used to estimate TH from anion, CL from TDS, anion from cation, CL from anion and vice versa and estimate CL from SAR. The results show that the concentration of calcium, potassium, sodium, chlorine, SAR, TH, PH is higher in the northwest and southeast, and as a result, the quality of water for various uses is lower, which is one of the possible reasons, in addition to the different formations of the earth. The geology and type of rocks (that is, the presence of gypsum materials and stones in the northwest and southeast parts or smaller amounts of these materials in the north, central and south parts (Owaisi, 2000)) can lead to higher temperature and drier environment which is the result of the higher rate of evaporation and transpiration as well as the type and density of land use (Laleh-Zari, 2019) [ABSTRACT FROM AUTHOR]