Salinity-based spatial evaluation of groundwater quality for agricultural use

Groundwater resources along with surface-water meet the needs of urban, industrial, and agricultural sectors, and therefore, their quality should be examined. Water quality indices are useful tools for aquifer management. In this research, the groundwater quality of the Miandoab Plain in the province of West Azerbaijan, northwest Iran, for agricultural purposes was investigated. For this purpose, the concentrations of the ions Mg2+, Ca2+, Na+, HCO3−, SO42−, Cl− and the pH level were measured. The indices effective salinity and potential salinity as well as sodium adsorption ratio and electrical conductivity were analyzed to evaluate the salinity. The geostatistical analysis was performed using the GS+ software, and the zoning maps of salinity hazard were prepared using ArcGIS. To prepare the maps, electrical conductivity (EC), effective salinity (ES), potential salinity (PS), and sodium adsorption ratio (SAR) as well as Mg2+, Ca2+, Na+, HCO3−, SO42, and Cl− were selected based on the semi-variogram values and cross-validation technique. The Cl− map was considered as the basis for preparing the groundwater quality maps of the region. The results showed that the groundwater quality in the east of the plain is suitable, in the central part can be recommended under constant supervision, and in the west is unsuitable for agriculture. In other words, according to the geography of the plain, the recharge area is the low-risk part of the plain and the salinity hazard increases toward the discharge area. The results can pave the way for the relevant organizations to plan for the agricultural and environmental sectors.


Introduction
Optimum use of groundwater requires proper management of exploitation which in turn needs scientific and principled knowledge of water resources in a region. Water quality indicators, as useful tools for water resources management, are usually intended for managers and decision-makers in a simple understandable way to transmit information on the quality and potential uses of a water body based on some criteria (Delgado et al. 2010). The Miandoab Plain is the main plain of Lake Urmia catchment and one of the agricultural hubs in the northwest of Iran. The turnover and economy of the plain is directly or indirectly dependent on the agricultural sector, most of the water needed for which is supplied from the groundwater resources. Almost half of the surface currents entering Lake Urmia pass through the Miandoab Plain, and 20% of the total groundwater abstraction in the catchment belongs to this plain. Therefore, proper management of the groundwater resources in the Miandoab Plain is necessary (Emami et al. 2018). Limited studies have been conducted on the groundwater of the Miandoab Plain: Emami et al. (2018) predicted the groundwater level of the Miandoab Plain using the artificial neural network (ANN), and the election (EA) and genetic (GA) algorithms. They concluded that the EA, with a root mean square error (RMSE) of 0.027, coefficient of correlation (R 2 ) of 0.93, and Nash-Sutcliffe efficiency coefficient (NSE) of 0.74, had a better performance in predicting the groundwater level than the ANN and GA. Noroozi-Qushbulaq et al. (2019) evaluated the vulnerability of the Miandoab Plain to nitrate using the GA. The optimized drastic (with the GA) had a higher correlation with nitrate, and therefore presented better results than the drastic model did. The optimized map showed that about 18,11,28,26, and 17 percent of the plain are located in the very low, low, medium, high, and very high vulnerable areas, respectively. Dehghanipour et al. (2019Dehghanipour et al. ( , 2020 presented and applied simulation-optimization approaches in the irrigated Miandoab Plain for identifying sustainable water management strategies. Results demonstrated the usefulness and flexibility of the methodology in identifying a range of potential water management strategies in complex irrigated endorheic basins like the Lake Urmia basin. Geostatistical method is one of the most accurate methods for estimating the values of variables at points where sampling has not been done. This method is based on the calculation of the weighted average of the values in the neighborhood of the estimation block. Geostatistical methods have widely been used in water engineering, some of which are as follows: Delgado et al. (2010), using the effective salinity and potential salinity, investigated the groundwater quality in Yucatan, Mexico. Gharbia et al. (2016) evaluated groundwater quality using GIS-based geostatistical algorithms. After analyzing the groundwater quality parameters, they created the maps of each parameter using the Kriging approach. Safiur Rahman et al. (2017) evaluated the water quality for sustainable agriculture in the district of Faridpur, Bangladesh. The irrigation water quality index showed that a large part of the samples was acceptable for irrigation. The spatial mapping of the water quality index showed that the groundwater was more suitable for irrigation in the north of Faridpur than in the central and southern parts. Taheri Tizro and Mohamadi (2019) used a geostatistical approach to evaluate the groundwater quality in the Zarin Abad Plain, Iran. They concluded that the electrical conductivity ranged between 480 and 6580 μS cm −1 , and the major cations and anions were ranked as Na + > Ca 2+ > Mg 2+ and SO 4 2-> Cl -> HCO 3 − . According to the values of RMSE and MAE, the Co-Kriging method was more optimum than the Kriging and inverse weighted distance strategies in evaluating the spatial variation of the groundwater quality parameters. Honarbakhsh et al. (2019) assessed the hydro-chemical quality of groundwater in semi-arid regions and mapped the qualitative parameters using ArcGIS. According to the Wilcox diagram, only 24% of the samples were in the C4-S4 class with high salinity and alkalinity. The maps showed that, due to the existence of dolomite and chalky formations in the south of the region, the groundwater in the north side was of better quality. ElKashouty (2019) studied the distribution of groundwater quality in the Nile Delta using geostatistical investigation (GIS) and concluded that the default Kriging was the best approach to map the hydro-geochemical parameters. Jeon et al. (2020) studied the trend of temporal variations in the groundwater quality in the Key sectors of Korea. For irrigation sustainability assessment, they used the US salinity laboratory (USSL), Wilcox diagram, electrical conductivity, sodium adsorption ration, sodium percent, and residual sodium carbonate. Based on the USSL and Wilcox diagram, about 98% and 83% of the samples were within the acceptable range for irrigation.
To our knowledge, despite the importance of quantity and quality of groundwater resources in the Miandoab Plain, there are a limited number of studies focusing on the distribution of groundwater quality in the plain from an agricultural perspective. On the other hand, no spatial analysis has yet been done on the groundwater resources in the plain, even though it can lead to better utilization of the resources and decrease in the soil destruction hazard. The aim of the present study is to analyze and assess the groundwater quality in the Miandoab Plain for agricultural uses based on geostatistical validation of major ions and salinity indices. This research was conducted in the Islamic Azad University of Tabriz in 2019-2020.

The study area
The alluvial aquifer of the Miandoab (Ghosha-Chay) Plain, as the fourth-most fertile plain in Iran, is the largest aquifer of Lake Urmia catchment. The plain is located inside the catchment in north-west Iran at longitude 36° 50′ to 37° 15′ east and latitude 45° 50′ to 46° 15′ north (Fig. 1). The aquifer is recharged from the north-east (the Mordaq-Chay River), east (the Lilan-Chay River), south-east (the Zarineh-Rood River), and south (the Simineh-Rood River) and discharged into the north-west (Lilan-Chay) and west (Zarineh-Rood and Simineh-Rood) of the plain. The Mordaq-Chay and Qoori-Chay rivers enter the plain from the north-east and east sides and immediately join the Lilan-Chay and Zarineh-Rood rivers. The Miandoab Plain has an area of 1256 km 2 and an average slope of 0.004. The plain is a part of the Alborz-Azerbaijan structural tectonic zone. The Zarineh irrigation and drainage network with an area of 586 km 2 is located in the plain. The Miandoab study area is divided into two main areas in terms of accessibility to the irrigation and drainage network: The areas located in the Miandoab network including Z 1 , Z 2 , Z 3 , Z 4 , Z 5 , Z 6 , Z 7 , Z 8 , Z 9-1 , Z 9-2 , Z 9-3 , Z 9-4 , Z 9-5 and the areas outside of the network including M 1 , M 2 , L 1 , L 2 , G, S, ZN, ZP. Therefore, the plain consists of the areas located in the Miandoab network and all the farmlands in L 2 , G, S, ZN and ZP (Jonubi et al. 2018).
The topography of the plain is uniform and smooth with only a few small hills in the central part. The plain is at an average elevation of 1292 m, with the lowest and highest points at the elevations of 1297 and 1400 m, respectively. The region has different geological formations. However, alluvial formations form the main part of the plain, which has emerged in the form of alluvial fans and fluvial deposits. The thickness of the alluvium varies across the plain, ranging from a few meters in the vicinity of the surrounding mountains to more than 150 m near the cities of Miandoab and Malekan. In the north of the plain, on the bed of the Mordaq-Chay River, the thickness of the alluvial deposits reaches up to 80 m. In the east, the thickness is about 30 m with thin layers of clay, but it even reaches up to 160 m due to the presence of the Lilan-Chay River (Emami et al. 2018).
Based on the Emberger's empirical climagram and using the data from the Miandoab synoptic station, the region has a cold semi-arid climate. As the aquifer is recharged by precipitation, the climate of the region is of great importance. The amount of precipitation in the region increases from the plain toward the heights. According to the thirty-year data set of 1988-2017, the average annual precipitation in the Miandoab study area is about 284 mm year −1 . The amount of potential evapotraspiration and pan evaporation are 742 and 1650 mm year −1 , respectively (Noroozi-Qushbulaq et al. 2019). According to the latest statistics, the amount of groundwater abstraction through 14,096 wells in the study area is equal to 350 MCM year −1 , 240 MCM of which belongs to the aquifer of the Miandoab Plain (Emami et al. 2018).

Method
To investigate the spatial and temporal variations of the groundwater quality parameters in the Miandoab Plain, the data set of 2002-2018 collected from 51 wells of the monitoring network by the Regional Water Company of West Azerbaijan was used. The data had been recorded in the form of laboratory results, solutes and necessary parameters including cations and anions. The emphasis was on the three cations sodium, calcium, and magnesium, and the four anions chlorine, sulfate, carbonate, and bicarbonate since these ions are more abundant than the other solutes in water. The long-term trend of the parameters are shown in Fig. 2.
First, all the data were normalized using the Kolmogorov-Smirnov test and then the classical statistical analysis of parameters was performed at a significance level of 5% in the SPSS22 software. If the data distribution did not follow the normal distribution, the logarithmic transformation was used. The Excel software was used to examine the presence or absence of a trend in the data. The empirical semi-variogram of the qualitative data was calculated using the GS + software and the best model was selected for each parameter. The interpolations were performed using the Kriging and inverse distance weighting methods. The interpolation methods were evaluated and validated using the cross-validation technique. To determine the best interpolation method for mapping the qualitative parameters, RMSE and ME as error criteria, and r as correlation coefficient were used; the method with the lowest RMSE and ME, and highest r was selected as the most appropriate one. Finally, using the ArcGIS software, the zoning maps of the groundwater salinity hazard were prepared.

Groundwater quality for agricultural use
Salinity is one of the most important indices in evaluating groundwater quality. One of the techniques for determining salinity is to measure electrical conductivity; the higher the amount of solutes in the water is, the higher the electrical conductivity is. The increase in salinity reduces the osmotic activity of plants and, as a result, the adsorption of water and nutrients from the soil (Saleh et al. 1999). Each crop, depending on its Fig. 1 The geographical location of the study area and sampling wells type, can tolerate salinity up to a threshold level above which the yields decrease per unit of increase in salinity (Maas and Hoffman 1977;Maas and Grattan 1999). In case of high salinity, some new indices such as effective salinity and potential salinity have been introduced (Delgado et al. 2010). Sodium adsorption ratio is also a suitable index to evaluate soil alkalinity hazard (Subramani et al. 2005;Kaur et al. 2017). There can also be chloride toxicity hazard in groundwater (Richards 1954). Hence, the electrical conductivity (EC), effective salinity (ES), potential salinity (PS), sodium adsorption ratio (SAR), and the ions affecting these indices, i.e. Mg 2+ , Ca 2+ , Na + , HCO 3 − , SO 4 2− and Cl − were chosen in the present study.
The Wilcox classification, also known as the US salinity laboratory method (Rafi Sharifabad et al. 2017), was used to determine the agricultural water quality. According to this classification and the values of EC and SAR, the agricultural water quality is classified into four categories: excellent, good, moderate, and not recommended (Table 1), and 16 classes (Table 2).
In the SAR classification ( Fig. 3), S 1 can be used for majority of soils which are less capable of reaching high exchangeable sodium content. S 2 can cause problems in fine-textured soils which do not contain gypsum and have a high cation exchange capacity, as well as when leaching is not possible. S 3 can cause the toxicity due to exchangeable sodium in irrigated soils; hence proper drainage, double leaching, and addition of organic matter would be necessary (however in soils containing large amounts of gypsum, there are less hazard). S 4 ; the soil is not suitable for irrigation, except in low or medium salinity conditions or when using gypsum as soil conditioner (Richards 1954;Delgado et al. 2010).
Electrical conductivity is an indirect measurement of the total concentration of salts dissolved in water, presented at the temperature of 25 °C.
In the EC classification (Fig. 3), C 1 can be used for a large number of crops and soils with low salinity hazard. C 2 can be applied to crops which are able to tolerate medium salinity. C 3 cannot be used for soils with poor drainage. Therefore, it is necessary to choose crops which are significantly tolerant of salinity, even when drainage is in good condition. C 4 is not usually suitable for irrigation, but if necessary, can be used for crops cultivated in highly permeable well-drained soils (Richards 1954;Palacios and Aceves 1970;Ayers and Wescott 1994;Delgado et al. 2010).
For the EC values greater than 250 µmhos cm −1 , ES and PS were also used since some salts can precipitate and in some others the relative concentration can increase due to evaporation (Palacios and Aceves 1970;Delgado et al. 2010).  Potential salinity is the estimation of risk of high salt concentrations due to the presence of Cl − and SO 4 2 , which can increase the osmotic potential of the soil solution when the available soil moisture content is less than 50%. The classification of water quality in terms of salinity based on the PS values is presented in Table 3. This index is defined as follows (Delgado et al. 2010). (2) Effective salinity is another index developed to consider relative solubility of various salts that are likely to exist in water. ES provides a more precise estimate of risk of osmotic pressure increase in soil solution when carbonate and bicarbonate concentrations are high (Kirda 1997). Under this condition, the calcium and magnesium carbonates and calcium sulfate precipitate, thus preventing an increase in the osmotic pressure in the solution. Water quality in terms of salinity based on ES is the same as the PS classification (Table 3). This index is described as follows (Delgado et al. 2010): where the sum of cations or the sum of anions means the use of the highest value (Palacios and Aceves 1970;Delgado et al. 2010).
Chlorides and some other ions, even in small amounts, can sometimes be toxic elements to plants and fruit trees. Water toxicity based on chloride, bicarbonate and sodium is presented in Table 4 (Ayers and Wescott 1994;Delgado et al. 2010).

Geostatistical analysis
Geostatistics is a subdiscipline of applied statistics that using information obtained from observed points provides Slightly saline, almost suitable for agriculture C 2 S 1 , C 2 S 2 , C 1 S 2 3 Saline, usable for agriculture after the necessary preparations C 3 S 3 , C 3 S 2 , C 3 S 1 , C 2 S 3 , C 1 S 3 4 Highly saline, hazardous to agriculture C 4 S 1 , C 4 S 2 , C 4 S 3 , C 4 S 4 , C 3 S 4 , C 2 S 4 , C 1 S 4

Sodium -Adsorption -Ratio (SAR )
Low -S1 Medium -S2 High -S3 Very High -S4 Low -C1 Medium -C2 High -C3 Very High -C4  a wide range of statistical estimators in order to estimate the desired characteristic at unobserved points. One of the tools for geostatistical studies is a function called variogram which allows the analysis of scale structure and intensity of spatial variation of regional variables. The GS + software is a comprehensive program which provides all geostatistical components including variogram analysis as well as mapping in an integrated and flexible way. Among the most important features of this software, which have also been used in the present study, are: spatial autocorrelation analysis, interpolation, and calculation of statistical parameters (Khosravi and Abbasi 2015).
.In geostatistics, to study the structure of the variability of the studied variable with respect to distance (spatial or temporal), it is necessary to establish an appropriate semi-variogram function. For this purpose, the sum of the squares of the difference between the point pairs that are at a known distance of h from each other should be calculated and drawn versus h (Sanches 2001 where Z(x i ) and Z(x i + h) represent the values of water quality at the point i and the other points isolated from x i at a distance of h, respectively.
Since the number of samples used in this study was 51, the isotropic semi-variogram was selected as the most suitable semi-variogram, because anisotropic semi-variograms require more samples. Among the common models, the one with a higher coefficient of determination (r 2 ), and lower nugget effect (C 0 ) and residual sum of squares (RSS) was selected as the best model to show the variogram.
The data were estimated using the Kriging and inverse distance weighting interpolation methods as Eqs. 8 and 9: where Z(x 0 ) is the optimum estimate of water quality, i is the optimum weight chosen to minimize the estimation variance of the sample i, Z(x i ) is the observed water quality, d i is the distance between x 0 and x i , p is the power of the parameter, and n is the number of observed data (Rafi Sharifabad et al. 2017).
The main difference between the above-mentioned interpolation methods is that the inverse distance weighting is a classical statistical method in which interpolation is based on the weight of points but Kriging is a geostatistical method expressing that in addition to the distance, the covariance of the data is also important.
There are different criteria for data validation, evaluation of the accuracy of the estimates and errors, and choosing the best interpolation method. In the present study, the cross-validation technique, mean error (ME), root mean square error (RMSE), and coefficient of correlation (r) were used to compare the observed and estimated values of the agricultural water quality (Hernandez-Stefanoni and Ponce-Hernandez 2006; Delgado et al. 2010).
The cross-validation technique is one of the most important and common techniques for data validation. In this method, one observation point is removed at each stage and then is estimated using the remaining observation points. This is repeated for all the observation points, so that at the end there will be an estimate for each observation point and finally, having the observed and estimated values, errors and deviations can be determined (Bagheri and Mohammadi 2009).
The mean error is used to determine the bias in the estimations according to Eq. 10 (Delgado et al. 2010): The root mean square error is a repeated measurement of the differences between the values estimated by a model and the observed values (Delgado et al. 2010). RMSE is a good criterion of accuracy and is used as an important parameter to show the precision of spatial analysis in ArcGIS (Nazarizadeh et al. 2006;Rafi Sharifabad et al. 2017). This criterion is calculated using Eq. 11 (Delgado et al. 2010): where n is the total number of samples, and Z x i and Z x i are, respectively, the estimated and observed values of water quality. Table 4 Water toxicity based on chloride, bicarbonate and sodium WR without restriction, S-M slight to moderate, S severe Water toxicity Cl − (mmol L −1 ) Hco 3 − (mmol L −1 ) Na + (mmol L −1 ) WR < 4 < 1.5 < 9.5 S-M 4-9 1.5-7.5 = 9.5 S > 9 > 7.5 > 9.5 The coefficient of correlation is a criterion of the relation between two variables and is denoted by r. The value of r varies between − 1 and + 1 (Delgado et al. 2010). This computation process is a comprehensive method for selecting the best ions and salinity indices in order to prepare zoning maps of spatial changes in groundwater quality for agricultural use (Delgado et al. 2010;Gharbia et al. 2016).

Chemical characterization of samples
To analyze the suitability of groundwater in the region for agriculture, the data set of 2002-2018 was classified using the Wilcox diagram. The results are presented in terms of frequency in Table 5. According to Table 5, saline water has the highest frequency, but fresh water and slightly saline water have no frequency during the period. Highly saline water has occurred in 37% of the cases.
The chloride concentration in 19 water samples (approximately 37% of the area of the plain) is less than 4 mmol L −1 , which has no restriction for agricultural use. The chloride concentration in 10 water samples (20% of the area) ranges from 4 to 9 mmol L −1 (low to moderate restriction). Finally, 22 samples (43% of the area) have a chloride concentration of more than 9 mmol L −1 (severe restriction).
The sodium content is less than 9.5 mmol L −1 (no restriction) in 25 water samples (49% of the total area) and more than 9.5 mmol L −1 (severe restriction) in 26 samples (51% of the area). However, no sodium percentage was found to fall into the low to moderate restriction category.
The bicarbonate concentration is 1.5-7.5 mmol L −1 (low to moderate restriction) in 40 samples (78% of the area) and more than 7.5 mmol L −1 (severe restriction) in 11 samples (22% of the area). Yet there is no bicarbonate percentage in the "no restriction" category.
The ions and salinity indices were analyzed using classical statistics and assessed (Table 6). According to Table 6, the standard deviation of parameters is about 1 indicate the data are close to average and have little scatter. The values of skewness and kurtosis were found to be in the ranges of − 0.18 to 0.68 and − 1.41 to 0.29, respectively, indicating the normal and symmetric distribution (Habibi 2016).

Spatial analysis
According to Tables 7 and 8, the nugget variance, structural variance, ME, RMSE, r 2 , and r were used to evaluate the estimated values.
The concentrations of all the ions using the Gaussian semivariogram with an r 2 of 0.5-0.98 are presented in Table 7 ( Fig. 4). The parameter [C/(C 0 + C)] × 100 indicates the ratio of structural variance to non-structural variance and determines what percentage of the data are structured; the more it tends to 1, the more structured the data become. The calculated structural variance is between 59.4 and 86.1%, suggesting a high spatial correlation of all the ions studied in the region. The values of nugget variance are small compared to those of structural variance, which also implies the identification of the spatial distribution pattern for these ions (Delgado et al. 2010). Highly saline, hazardous to agriculture 37% The indices ES, PS, and EC were explained by the Gaussian model with an r 2 of 0.98, and SAR was explained by the linear to sill model with an r 2 of 0.94 (Fig. 5). The values of structural variance and nugget variance are 82.7-89.8% and 18-53.1%, respectively (Table 7). Accordingly, the spatial autocorrelation in the agricultural water quality indices was clearly identified by determining the Gaussian and linear models appropriate to the empirical semi-variograms (Fig. 5).
The cross-validation of the interpolations for the ions showed that r, ME, and RMSE of the observed and estimated values range between 0.22 and 0.71, 0.18 and 0.66, and 0.22 and 0.84, respectively. For EC, PS, ES, and SAR, these criteria are in the range of 0.62 to 0.67, 0.43 to 0.65, and 0.54 to 0.81, respectively (Table 8).
In all the parameters, the ME and RMSE values indicated that the errors in the interpolation models are small, and therefore acceptable. The correlation coefficient of Cl − is high and HCO 3 − , ES, SAR, Na + , PS, and EC have good correlation coefficients, respectively, so were selected as criteria for preparing the groundwater salinity zoning maps of the Miandoab Plain. However, no map was prepared for Ca 2+ , Mg 2+ , and SO 4 2− as they have correlation coefficients of less than 0.49 (Figs. 6 and 7).

The spatial distribution of the groundwater quality for agricultural use
The chloride map, due to having the best indicators of spatial analysis, was selected to identify the areas with different groundwater qualities. This map has three zones with no restriction, low to moderate restriction, and severe restriction for agricultural use, according to which the Miandoab Plain was classified into three zones: the east, center, and west of the plain.

Zone (1); east plain
According to the zoning maps (Figs. 6 and 7) of the east side of the plain, it was found that the EC and ES values are indicators of moderate water quality, the PS values with two ranges represent good and moderate quality, and SAR shows  (Table 9). In this zone, the concentrations of sodium and bicarbonate are low, indicating no restriction and low to moderate restriction, respectively, in terms of water toxicity. The chloride content in both ranges indicates low concentration, and thus low toxicity, representing no restriction and low to moderate restriction (Table 9). Therefore, the water quality in the east side of the plain is in good to moderate condition and is recommended for agricultural use.

Zone (2); central plain
The zoning maps (Figs. 6 and 7) of the central parts of the plain show moderate and unrecommendable water quality based on EC and ES values with two ranges. The SAR values indicate both excellent and good quality, and the PS values represent moderate quality. The sodium concentration is in two ranges with no restriction and severe restriction. The chloride concentration is almost low, indicating low to moderate toxicity. The bicarbonate concentration is also in two ranges with toxicity of higher than 3.5 mmol L −1 , representing moderate to severe restriction. Overall, in order to use the groundwater for agricultural purposes in the central parts of the plain, in addition to monitoring the soil salinity, the aquifer in the east and west of the plain should be taken into account and the quality indices and other necessary measures should be applied (Table 9).

Zone (3); west plain
The zoning maps (Figs. 6 and 7) of the west parts of the plain demonstrate low water quality based on the values of EC, ES, and PS (Table 9). Owing to the high values of these three parameters, the water quality cannot be recommended (Richards 1954;Palacios and Aceves 1970). In this zone, the concentrations of sodium and chloride are very high, resulting Fig. 4 The empirical semivariograms for Calcium, Magnesium, Sodium, Chloride, Sulfate and Bicarbonate in the groundwater in severe toxicity. The bicarbonate content with two ranges has a toxicity level of higher than 3.5 mmol L −1 , indicating moderate and severe restriction. However, in the west of the plain, only in terms of SAR the water quality is good and excellent with no restriction (Table 9). In general, the results show that in this zone, due to the high salinity and the high concentrations of chloride, bicarbonate, and sodium, the toxicity hazard to crops is severe. Therefore, the water in the west plain is of low quality and is not recommended for agricultural use. This can be a warning sign of the extreme exploitation of groundwater. According to the zoning maps, chloride shows low to moderate toxicity in the aquifer recharge area, which considering the parameters of potential salinity (Eq. 2), this can also be the reason for the moderate potential salinity in this area. Similarly, toward the discharge area, with increasing the chloride toxicity, the potential salinity also increases (Figs. 6a and 7d). Since bicarbonate is involved in the calculation of effective salinity (Eqs. 3-6), the bicarbonate toxicity at the entrance compared to the end of the plain can indicate effective salinity in the recharge and drainage areas to some extent (Figs. 6b and 7c).
According to the general trend of changes in groundwater, i.e. the Chebotarev sequence (1955), groundwater at the entrance (recharge area) is bicarbonate, which gradually changes to sulfate in the middle of the plain and finally, to chloride at the end (discharge area). Thus, considering the concepts of effective salinity and potential salinity, the bicarbonate and chloride contents Fig. 5 The empirical semivariograms for the Electrical Conductivity, Sodium Adsorption Ratio, Effective Salinity and Potential Salinity of the groundwater Fig. 6 The zoning maps of the ions contents in the groundwater respectively in the recharge and discharge areas (Fig. 6a, b), confirm the changes in the effective salinity and potential salinity (Fig. 7c, d).

Conclusion
The spatial and temporal changes in the groundwater quality in the Miandoab Plain were investigated using the data set of 2002-2018 from 51 wells. To evaluate the quality of groundwater for agricultural use, two new indices, effective salinity and potential salinity, along with EC and SAR and Wilcox diagram were analyzed. The classical and geostatistical analysis were performed using SPSS and GS + , respectively. The zoning maps of the groundwater salinity hazard were prepared using ArcGIS. According to the cross-validation technique, the chloride semi-variogram obtained by the Gaussian model displayed the best spatial distribution of EC, ES, PS, and SAR as well as HCO 3 − and Na + in the plain. The zoning maps of EC, ES, PS, SAR, Cl − , HCO 3 − , and Na + showed that the Miandoab Plain could be divided three zones with no restriction (zone (1); east plain), low to moderate restriction (zone Fig. 7 The zoning maps of electrical conductivity, sodium adsorption ratio, effective salinity and potential salinity in the groundwater Table 9 The zones of the groundwater quality for agricultural use based on the hazards of toxicity, salinity and sodicity (2); central plain), and severe restriction for agricultural use (zone (3); west plain). In other words, groundwater salinity and toxicity increase from the recharge areas (the north-east, east, south-east, and south) toward the discharge areas (the north-west, west, and partly south-west) of the plain. It can be concluded that salinity is higher in the west of the plain, which can partly be due to the existence of Lake Urmia in the west side. Since the aquifer is discharged into Lake Urmia, which is in turn the cause of the salinity, with the increase of groundwater abstraction from the plain, the way for the saline water to flow from the lake to the aquifer is facilitated, which in turn, intensifies the salinity. Therefore, the recharge area (east) is the low-risk side, and the discharge area (west) is the high-risk side of the plain. This can be a warning sign of the extreme exploitation of groundwater and should be considered in the management of water resources.
Funding This research received no external funding.

Declarations
Conflict of interest Authors certify that there is no conflict of interest in this research.
Ethical approval This article does not contain any studies with human participants or animals performed by any of the authors.