Introduction

Quantity and quality of groundwater are worsening due to the increase in urbanism and human activities (Shirazi et al. 2013; Alam et al. 2014; Duarte et al. 2015; Shekhar et al. 2015). In Torqabeh and Shandiz County, utilization of groundwater resources is highly important due to the county’s small area, low annual precipitation, and shortage of water resources. Considering the limited access to water resources and water quality problems, it is necessary to use quick techniques for identification and evaluation of groundwater resource vulnerability in this area. Vulnerability maps consider the geology of surface and also under the ground lands and it is a possibility of infiltration and spread of pollutants from the ground surface into the groundwater system (Piscopo 2001). Groundwater vulnerability, e.g., DRASTIC method, considers hydro-geological parameters (hydrological, hydrogeological, and geological factors). Many methods have been proposed for vulnerability assessment, including process-based methods, statistical methods, and overlay methods. These methods were used in research by Rahman (2008), Magesh and Chandrasekar (2013), Al-Shatnawi et al. (2014), Mogaji et al. (2014), Alaya et al. (2014), Sadat-Noori et al. (2014), Davraz and Özdemir (2014), Chandrasekar et al. (2014), Narmada et al. (2015), Jhariya et al. (2019), and Nadiri et al. (2019). To draw groundwater vulnerability maps, the DRASTIC empirical method is often used in the assessment of groundwater vulnerability (Shirazi et al. 2012, 2013; Jilali et al. 2015; Shekhar et al. 2015; Ayed et al. 2017). The DRASTICindex was designed by Aller et al. (1987) in collaboration with the US Environmental Protection Agency and the US National Water Association, only with consider the hydrological parameters to assess the vulnerability of groundwater to pollutants. Then, this index was developed by many researchers with considering geological, hydrological, weather, and other parameters for assessment of groundwater vulnerability (Chandoul et al. 2014). This method assumes that some major factors involved in controlling groundwater vulnerability can be weighted. Since the data required by the DRASTIC method are easily accessible, this method is extensively used in many countries for groundwater vulnerability assessment (Babiker et al. 2005; Chitsazan and Akhtari 2009; Huan et al. 2012; Yin et al. 2013; Kaliraj et al. 2014; Neshat et al. 2014; Duarte et al. 2015; Jilali et al. 2015; Singaraja 2015; Su et al. 2015; Ghosh et al. 2015; Ayed et al. 2017). The geographic information system (GIS) can be used to obtain more reliable results more quickly in determining the value of the DRASTICindex for groundwater vulnerability assessment (Ahmed 2009; Al-Adamat et al. 2003; Baalousha 2011; Albuquerque et al. 2013; Magesh and Chandrasekar 2013; Chandrasekar et al. 2014; Davraz and Özdemir 2014; Kaliraj et al. 2014; Varol and Davraz 2014; Zhang et al. 2014; Mondal et al. 2019). Kura et al. (2015) used the GIS system and the spatial relationship of groundwater quality with geological and topographic parameters, land coverage, and sources of pollution for assessing groundwater pollution. Some of the studies that employed the DRASTIC method to assess groundwater vulnerability include the research by Saidi et al. (2010), Uhan et al. (2011), Mishima et al. (2011), Huan et al. (2012), Yin et al. (2013), Gupta (2014), Mogaji et al. (2014), Neshat et al. (2014), Varol and Davraz (2014), Wu et al. (2014), Narmada et al. (2015), Shekhar et al. (2015), Su et al. (2015), Ayed et al. (2017), and Singha et al. (2019).

This research was conducted in three phases: development of spatial databases, analysis of spatial data, and data merging. The spatial databases of the DRASTIC layers (depth of water table, net recharge, aquifer media, soil media, topography, the impact of the vadose zone, and hydraulic conductivity) were obtained from different resources. After digitizing, the related maps were prepared, errors were corrected, and descriptive information was added to the maps. Finally, the maps were overlaid in GIS in a spatial format.

The objectives of the paper are as follows: (1) to assess vulnerable zones modified DRASTICindex by applying a novel methodology, including DRASTIC—land use in Torghabeh-Shandiz of northeast Iran. Nitrate is selected as the main contaminant introduced into the groundwater environment of the study area by human activities, and as a representative indicator of groundwater quality degradation; (2) to validate the accuracy of the modified vulnerability index; and (3) was adopted the analysis of variance of F statistics to verify the groundwater vulnerability results.

Study area

Torqabeh and Shandiz County are located in the Razavi Khorasan Province with an area of 1184.82 km2. It is situated on the 59° 0.0′ to 59° 40′ eastern longitudes and 36° 0.0′ to 36° 30′ northern latitude (Fig. 1). The climate of the region is cold, arid, and naturally mountainous and plain areas. Different parts of the region have high erosion potential due to high slope and low soil depth and low-density vegetation with low rainfall. The vegetation of the area is located in the Iranian-Turanian vegetation zone based on the divisions of the Iranian Geographic Map. The area has a variety of lowlands and high plateaus, and the climatic conditions are almost similar, with low rainfall and a long dry season and high thermal fluctuation. The average temperature and annual precipitation for a 25-year period (1993 to 2018 years) are 13.3 °C and 234 mm, respectively (Irimo 2018).

Fig. 1
figure 1

Sampling points in the study area of Torghabeh-Shandiz of Khorasan County, Iran

Full size image

Materials and methods

DRASTICindex

DRASTICindex is developed for the assessment of groundwater vulnerability on the regional scale. The variables involved in this index are classified into three main categories: land surface factors, unsaturated zone factors, and vadose zone. This index is based on the following four hypotheses (Adnan et al. 2018): (1) ground surface pollution, (2) percolation of pollution into groundwater through precipitation, (3) water causes the movement of pollution, (4) the minimum area of the DRASTIC-assessed zone is 0.4 km2. The DRASTIC index’s parameters include the following: depth of water table (D), net recharge (R), aquifer media (A), soil media (S), topography (T), the impact of vadose zone (I), and hydraulic conductivity (C) (Pacheco et al. 2015). To assess the degree of groundwater pollution in the region using the DRASTICindex, the Spatial Analyst tool extension in ArcGIS 10.3 was used. The DRASTICindex is calculated using Eq. (1).

$$ {\mathrm{DRASTIC}}_{\mathrm{index}}={D}_{\mathrm{r}}{D}_{\mathrm{w}}+{R}_{\mathrm{r}}{R}_{\mathrm{w}}+{A}_{\mathrm{r}}{A}_{\mathrm{w}}+{S}_{\mathrm{r}}{S}_{\mathrm{w}}+{T}_{\mathrm{r}}{T}_{\mathrm{w}}+{I}_{\mathrm{r}}{I}_{\mathrm{w}}+{C}_{\mathrm{r}}{C}_{\mathrm{w}} $$
(1)

where (D) depth to water level, (R) recharge, (A) aquifer media, (S) soil media, (T) topography, (I) impact of the vadose zone, (C) hydraulic conductivity, and the subscripts r and w are the corresponding ratings and weights, respectively (Jhariya et al. 2019; Nadiri et al. 2019).

In the DRASTICindex, the amount of pollutants that reach to the groundwater was determined, and higher values of this index indicate a higher degree of pollution (ranked 1 the lowest and 10 the highest vulnerability of groundwater pollution). Afterwards, a weight between 1 and 5 was assigned to each of the seven parameters based on their relative qualitative importance for the potential for the conveyance of pollution in the groundwater system (Table 1). Finally, by multiplying the weight and rank of each of the seven layers influencing pollution in ArcGIS, the DRASTICindex was obtained an overlay using Eq. (1). The higher DRASTICindex shows more extensive pollution of groundwater resources in the area. The methodology chart for groundwater vulnerability study is shown in Fig. 2.

Table 1 Parameters of the DRASTICindex by Pacheco et al. (2015)
Full size table
Fig. 2
figure 2

Flowchart of the method for groundwater vulnerability assessment

Full size image

Results and discussion

Spatial Analyst tools of GIS environment have been used to determine the DRASTIC vulnerability index by incorporating all the thematic layers and reclassifying by assigning appropriate ratings and weights.

Depth-to-water table

There is an inverse relationship between the depth to the water table and pollution of groundwaters. That is to say, with an increase in depth to the water table, pollutants reach the water table more slowly and the possibility of diffusion/dilution or absorption of pollutants in soil grows (Shekhar et al. 2015). The depth-to-water table parameter was derived from water level data taken from 35 piezometric wells located in the study area. The depth-to-water table from ground level point information was interpolated, using the inverse distance weighting (IDW) method in ArcGIS 10.3, to derive the depth to the groundwater table surface. Thus, shallower the water depth, the more vulnerable to aquifer pollution. Afterwards, groundwater depth was ranked (Aller et al. 1987) by assigning a rank between 1 and 10 to each depth based on their pollution potentials (1 refers to the lowest and 10 refers to the highest vulnerability for groundwater pollution). For instance, ranks 10 and 9 were assigned to depths smaller than 20 m (< 20 m) and depths between 20 and 40 m (20–40 m), respectively. Other depths were ranked based on Table 2, while a weight of 5 was assumed for this layer. Figure 3 a shows the result of the depth to water map.

Table 2 The ranges, ratings, and weight of the DRASTIC quantitative parameters (source: Aller et al. 1987; Venkatesan et al. 2019)
Full size table
Fig. 3
figure 3figure 3

Depth of water level (a), soil permeability (b), net recharge (c), aquifer media (d), soil media (e), topography (f), impact of the vadose zone (g), and hydraulic conductivity (h)

Full size image

Net recharge

Net recharge is the total amount of water per unit area, in mm per year, which reaches the water table. Therefore, the recharge is an essential vehicle for percolating and transporting pollutants in the vadose zone to the saturated zone (Shirazi et al. 2013). When the recharge is more, the higher the chance of pollutants to reach the water table. Since pollutants move further into the groundwater with increased mobility of water, the net recharge is a crucial factor. The amount of recharge that percolates into the aquifer depends on the path of flow from the land surface (Kura et al. 2015). The net recharge was calculated by using the Piscopo (2001) formula shown in Eq. 2.

$$ R=\mathrm{In}+R+S $$
(2)

where R = Recharge layer, In = Infiltrations (soil permeability index), R = Rainfall index (mm), and S = Slope index.

The slope index was determined using the digital elevation model (DEM), provided by the National Cartographic Center of Iran, was generated with a spatial resolution of 30 m. After deriving the slopes, the slope was reclassified according to the criteria given in Table 2. The influence map of the soil was obtained from the Environment Organization and was digitized in ArcGIS 10.3. In order to prepare the rainfall layer, a 25-year period (1993 to 2018 years) data of the Golmakan Meteorological Station were used based on the relationship between the amount of rainfall and the height of the digital elevation model (Fig. 3b). The slope, precipitation, and soil type parameters were classified using the criteria in Table 2. Finally, the aforementioned three layers were overlaid in ArcGIS and the net recharge layer was obtained. Based on Table 2, a weight of 4 was assumed for the net recharge layer. The influence map and net recharge are shown in Fig. 3c.

Aquifer media

Aquifer refers to the portion of ground capable to yield water in pores. It is the saturated zone material, which controls the polluted attenuation processes which determine the flow rate. The larger the grain size, the more fractures or openings within the aquifer, and the higher permeability (Venkatesan et al. 2019). To prepare the aquifer media, the drilling log information of piezometric and harvest wells in the region was used. After collecting the wells’ log information, the vertical lithological column information on several deep wells was determined, and the dominant type in the phreatic zone and aquifer media was obtained in ArcGIS 10.3 using IDW method. Afterwards, the result was ranked based on the aquifer media formation type. The weight of this layer was calculated based on Table 2. The weighted map of the aquifer media is shown in Fig. 3d.

Soil media

Soil significantly affects the amount of recharge penetrating into the water table and flow of pollutant. Soil pollution potential depends on soil type/material, particle size, and soil potential for pollution absorption. With a decrease in the diameter of soil particles, soil permeability decreases and soil thickness escalates. As a result, groundwater pollution potential declines (Al-Shatnawi et al. 2014). Using the drilling log information piezometric and harvest wells in the region, the map of soil media and the material was prepared to a depth of 2 m in ArcGIS 10.3 and the layer was weighted based on soil material (Table 2). The map of this layer is shown in Fig. 3e.

Topography

Topography has been specified as an indication of whether a pollutant has run off or remains on the surface extended enough to infiltrate into the groundwater. The slope can be defined as the measure of the steepness of an area. The area has a low slope leads to storage of water that gives more time for water to percolate as compared with areas having a large slope. The areas of the low slope have more vulnerability to groundwater pollution (Fernandes et al. 2014). The slope was determined from DEM of 30 m by using ArcGIS 10.3. Finally, based on Table 2, the rank of slope levels and weight of the topography layer were calculated. Figure 3f shows the topographic level map (slope map) of the region. The slope level map indicates that the region’s slope varies from 0 to more than 35%. The slope percentage of the region was classified into the following 5 classes: 0–5%, 5–15%, 15–25%, 25–35%, and > 35%. The slope of most of the area’s lands was in the 5–15% class (457 km2), and groundwaters in this slope class were more prone to pollution.

Vadose zone

The vadose zone starts from the surface soil zone and continues to the water table. This parameter shows the influence of the vadose zone over the water surface and controls the path of flow of pollutants into the aquifer. The vadose zone plays a significant role in reducing groundwater pollution due to processes such as biological decomposition, treatment, mechanical properties, and chemical reactions on the way of pollutants of this layer (Shirazi et al. 2013). Type of the vadose zone determines recharge properties and dilution of materials below soil profile and above the water table. This layer depends on soil material. In the DRASTICindex, it is assumed that the settings of the vadose zone largely affect pollutants because pollutants find the chance to be absorbed and diluted in this region before reaching the water table (Nazzal et al. 2019). Area of levels, rank, and weight of the vadose zone layer was calculated based on Table 2. Figure 3g shows the map of the vadose zone layer’s levels. The groundwater resources in the west and center of the region are more vulnerable to pollution. Therefore, the rank and weight assigned to these areas are large (Table 2).

Hydraulic conductivity

This layer deals with water table influence or potential of water table contents for transferring water or soluble materials. Hydraulic conductivity is the degree of the flow of groundwater under the media’s hydraulic slope (Shirazi et al. 2013). With an increase in hydraulic conductivity, the possibility of a flow of pollutants increases, and hydraulic conductivity depends on soil texture. As the size of soil grains increases, hydraulic conductivity escalates. Hydraulic conductivity is the factor controlling movement and retention time of pollutants from the moment they penetrate into soil surface at the moment they travel inside the water table. Hence, a rise in hydraulic conductivity increases pollution potential. Since this layer depends on soil texture, its rank was determined based on the soil map (Shekhar et al. 2015). The area and weight of the hydraulic conductivity layer’s levels were determined by Table 2, while the map of this layer is depicted in Fig. 3h. Hydraulic conductivity is higher in the west and center of the study area. Hence, these areas are more vulnerable and the soil in this region lacks vegetation and coarse-grained texture.

Vulnerability map with DRASTICindex

The seven resultant parameter map layers were overlaid using the weighted sum method of multiplying each parameter of weight with their corresponding rate is shown in Table 2, summing them to get the final vulnerability index. The study area was divided into five vulnerability classes using the Natural break algorithm from power tool of ArcGIS, which are ranging between a minimum value of 140 and a maximum value of 200. These classes are very low, low, moderate vulnerable, highly vulnerable, and very high vulnerable, according to the DRASTICindex and vulnerability class. The vulnerability zone map is shown in Fig. 4. DRASTIC vulnerable map shows that about 24.1 and 15.3% of the area status between high vulnerability to very high vulnerability of groundwater pollution and the remaining 24.2, 27.7, and 8.8% area are occupied by moderate vulnerability to the low and very low vulnerability of the pollution zone (Table 3). The high-vulnerability zones are present in the northeast and east parts of the watershed. This may be due to higher net recharge (12–15 cm/year) in these regions which carries the pollution of the surface along with the percolating water to the deeper aquifer zones which induces high pollution in these regions. The soil type also helps in easy infiltration of the contaminants to reach the groundwater of this area. The slope of the area also accelerates the pollution process as the pollution gets more time to infiltrate the flat topography. The hydraulic conductivity in the northeast and east region also increases the chances of groundwater vulnerability to pollution in this region. The very high vulnerability zone covers an area of 181.8 km2 and is present towards the northeast and east regions.

Fig. 4
figure 4

Groundwater vulnerability map with the DRASTICindex

Full size image
Table 3 Area under vulnerability (DRASTICindex) to groundwater pollution
Full size table

Land use extraction

Remote sensing image with the clear-sky condition was downloaded (from https://landsat.usgs.gov). The Landsat-8 Operational Land Images acquired in 2018 was used in this research. The remote sensing (RS) data were processed with the following steps. First, a geometric correction was made for the RS data using the polynomial method based on ground control points. To unify the data resolution, all the data were resampled to 30 × 30 m. After radiation correction and atmospheric correction, the False Color Composite (FCC) was used for diagnosing feature. These operations were conducted on ERDAS Imagine 9.3 Software (Rawat and Kumar 2015). Then, seven typical categories of land use: gardens, residential area, irrigated agriculture, dry farming, wetland, forests, woodlands and rangelands (low, semi-dense, and dense rangeland) were extracted using three spectral index-based methods. The cell size of the extracted land use was set at 30 × 30 m. The water bodies inside the study area were extracted using the modified normalized difference water index (MNDWI) (Guo et al. 2019):

$$ \mathrm{MNDWI}=\frac{\mathrm{Green}-{\mathrm{MIR}}_1}{\mathrm{Green}+{\mathrm{MIR}}_1} $$
(3)

where Green and MIR1 refer to the green band and a middle infrared band of the RS image, respectively. Next, a manual threshold for separating all of the bodies of water from other areas was determined. This method was modified based on the MNDWI method. The results contained more information and eliminated the terrain influence (Acharya et al. 2018). The area covered with vegetation was extracted using the normalized difference vegetation index (NDVI) method after masking out water areas (Huete et al. 2002):

$$ \mathrm{NDVI}=\frac{\left(\mathrm{NIR}-\mathrm{RED}\right)}{\left(\mathrm{NIR}+\mathrm{RED}\right)} $$
(4)

where NIR and RED refer to a near-infrared band and a red band of the RS image, respectively. The threshold (0.255) was determined after a thorough visual interpretation of the study area. To work out the land use classification, supervised classification method with maximum likelihood algorithm was applied in the ERDAS Imagine 9.3 Software (Fig. 5).

Fig. 5
figure 5

Land use maps

Full size image

Land use and human activities, a considerable influence on the groundwater vulnerability of most of the area. A region’s pollution potential varies by its land use pattern. Hydrologic parameters are influenced by land uses. Agricultural pesticide, drilling well, harvest from mines, landfills, and commercial and industrial wastes may change hydrologic parameters. Land use classifications of the region (Table 4) show that a major use in the region is rangelands. Agricultural lands are the second most extensive lands in this region, while other land uses include industrial and urban land uses. Groundwater in agricultural areas is more vulnerable to nitrate concentration (Mishima et al. 2011; Uhan et al. 2011). Nitrate distribution in groundwater mainly depends on soil dynamics, recharge level, the flow of groundwater, and nitrogen load on land (Uhan et al. 2011). The presence of nitrate in groundwater is influenced by agricultural and human activities and affects the region’s pollution degree.

Table 4 Types, area, rank, and weight of land use
Full size table

Vulnerability map with modified DRASTICindex

To assess the vulnerability of groundwater, the DRASTIC map was combined with the land use map. In this research, the vulnerability map was prepared by adding the land use parameter and combining it with the DRASTICindex. This combination is known as the modified DRASTICindex. The effect of agricultural, industrial, and urban activities on the groundwater vulnerability is greatly conformed to the vulnerability map. To prepare the vulnerability map, the land use map was ranked and weighted based on Table 4 (Saidi et al. 2010). Afterwards, the land use map was converted to the raster format and was multiplied by the weight of the parameter (Lw = 5). To create a spatial relationship between land use and the DRASTIC map, the land use map and DRASTIC map were overlaid. The final outcome of the overlaying of the DRASTICindex (DI) with land use was the modified DRASTICindex (MDI) (Table 5), which was calculated using Eq. (5).

Table 5 Area under vulnerability (modified DRASTICindex) to groundwater pollution
Full size table
$$ {\mathrm{MD}}_{\mathrm{index}}={D}_{\mathrm{index}}+{L}_{\mathrm{r}}{L}_{\mathrm{w}} $$
(5)

where Lr and Lw denote the rank and weight of the land use parameter, respectively. In addition, Dindex is the DRASTICindex and MDindex is the modified DRASTICindex.

The vulnerability map suggests that parts of the study area and different types of human activities have the highest contribution to the groundwater vulnerability. In the present research, Gupta (2014) and Kura et al. (2015) classification was used in the modified DRASTICindex to classify the vulnerability of the Torqabeh-Shandiz aquifer. Moreover, a combination of the DRASTICindex and land uses was used in this study to assess the groundwater vulnerability. The vulnerability map was divided into five categories: very low (< 100), low (100–139), moderate (140–175), high (175–200), and very high (200 <) vulnerability (Table 5). Figure 6 shows the modified DRASTIC aquifer vulnerability map. The results of the analysis indicated that 25.9% and 20.8% of the areas were in the high and very high vulnerability classes, 44% moderate vulnerability, and 9.3% low and very low vulnerability. A comparison of the DRASTIC map with the land use map also indicated that the quality of groundwater is immensely influenced by agricultural, industrial, and urban activities.

Fig. 6
figure 6

Modified DRASTIC aquifer vulnerability map

Full size image

The result of the difference area% between DRASTIC and modified DRASTIC methods shows that the percentage of low vulnerability areas have decreased by 27.4%, by applying this modified method. In contrast, the areas of moderate to very high vulnerability increased by 19.8% and 5.4%, respectively, this increase is caused by human activities (agricultural, industrial, and urban activities) (Table 6).

Table 6 Difference area% between DRASTIC and modified DRASTICindexes
Full size table

The results of the DRASTIC and modified DRASTICindex showed that a large part of the groundwater was in the moderate to high vulnerability classes. The results of the DRASTICindex indicated that 8.8%, 27.7%, 24.2%, 24.1%, and 15.3% of the area were in the very low, low, moderate, high, and very high vulnerability classes, respectively. The results of the modified DRASTICindex indicated that 9.0%, 0.3%, 44.0%, 25.9%, and 20.8% of the area were in the very low, low, moderate, high, and very high vulnerability classes, respectively. In view of the results of DRASTIC indexing, the vulnerability of the Torqabeh-Shandiz aquifer was mainly in the low to high range in both indexes. Pollution level and pollution potential of the central part of the aquifer were in the low vulnerability class, whereas the northeast and east of the aquifer were in the high to very high classes.

Validation of groundwater vulnerability to nitrate mapping by MDindex

To validate the vulnerability map prepared with the modified DRASTICindex in the field conditions, a sample of groundwater was obtained to analyze nitrate (NO3−) is selected to be the typical pollutant not only because it constitutes the main contaminant introduced into the groundwater environment of the study area by human activities, but also because it has been proposed as a representative indicator of groundwater quality degradation. Using a DR/4000 Hach spectrophotometer, the nitrate levels of the sample were measured (Huan et al. 2012). The nitrate concentration map was prepared in ArcGIS 10.3 using the IDW method (Fig. 7). To evaluate the accuracy of the nitrate map, nitrate concentration in 26 wells was sampled, and the correlation between actual (field measurement) and estimated values in sample points was determined (Fig. 8). Thus, the vulnerability map obtained using the modified DRASTICindex can be used to enforce land management development policies (Baalousha 2011; Tiwari et al. 2016). By adding nitrate concentration matched with the modified DRASTIC map, it was found that all of the points with high nitrate concentration were in the very polluted class, and this finding can approve precision and accuracy of the method. The pollutions are imparted to the groundwater due to the heavy application of fertilizer for agriculture.

Fig. 7
figure 7

Nitrate concentrations map using the IDW interpolation method

Full size image
Fig. 8
figure 8

Relationship of measured and estimated value of nitrate concentration

Full size image

The corresponding relationship between the nitrate concentration and the vulnerability classes, analysis of variance (ANOVA) of F statistic was adopted to verify the groundwater vulnerability results (McLay et al. 2001). The ANOVA of F statistic is, the large overlap there is between the nitrate values in different vulnerability classes. Formally, the F statistic is based on the significance test of the regression coefficient to judge the fitting degree of multivariate linear regression. Variance (ANOVA) F statistic can be calculated by Eq. (6).

$$ F=\frac{\mathrm{MST}}{\mathrm{MSE}}=\frac{\mathrm{SST}/k-1}{\mathrm{SSE}/n-\mathrm{k}} $$
(6)

where MST and MSE are respectively the mean square for treatment and mean square for error. SST and SSE are respectively the sum of squares for treatment and the sum of squares for error. k-1 and n-k are respectively the freedom degree for treatment and freedom degree for error (Huan et al. 2012).

The validation of the initial application of the MDindex depended on the correlation between the vulnerability index and nitrate concentration values of 51 groundwater samples collected in the common water period of 2018. The maximum, minimum, and average concentrations of NO3− were respectively 85 mg/L, 1.2 mg/L, and 23.32 mg/L (Fig. 9). A summary of the descriptive statistics of the distribution of nitrate concentrations in the study area is presented in Table 7. The result of the distribution of nitrate values in different vulnerability classes was not normal (K-S < 0.05 data is not normal). Then, the correlation between the MDindex and nitrate concentration values was expressed by the Spearman rank correlation factor which was high, in the order of 0.964 (Table 8).

Fig. 9
figure 9

Map of nitrate concentration matched with the concentration of pollution of groundwaters in the Modified DRASTICindex

Full size image
Table 7 Descriptive statistics of nitrate data in sampling points
Full size table
Table 8 Correlation factors between nitrate concentration and modified DRASTICindex
Full size table

By comparing the analysis of variance (ANOVA) F statistic with the Duncan test in Table 9, it can be noticed that the nitrate concentration classifications and the vulnerability classifications correspond to each other well. These observations led to the conclusion that the modified DRASTICindex reflect the actual groundwater vulnerability.

Table 9 ANOVA statistics by Duncan test for comparing modified DRASTICindex classes and nitrate concentration
Full size table

Conclusion

In this research, the vulnerability zoning map of the Torqabeh-Shandiz aquifer was prepared using the DRASTIC and modified DRASTICindex. In the DRASTICindex, seven parameters (depth to water level, net recharge, aquifer media, soil media, topography, the impact of vadose zone, and hydraulic conductivity) were used. By assigning weights and ranks to these seven factors, the layers were combined with GIS. In the modified index, in addition to the assignment of weights and ranks to the seven factors, the vulnerability map of the aquifer was obtained using the weights and ranks of land use classes and combined with the DRASTICindex in GIS.

A comparison of areas with high and very high vulnerability in the DRASTIC and modified DRASTICindex map indicated that these areas increased by 7.2%. This indicates that the modified DRASTICindex more attention to high and very high vulnerability classes. A comparison of the modified DRASTIC map to the land use map and nitrate concentration indicates that the quality of the groundwater is largely influenced by agricultural activities. Furthermore, the relatively low and high vulnerability classes corresponded respectively to the lowest and highest NO3− concentrations. The average concentrations of NO3− across vulnerability classes increased from low vulnerability to high vulnerability. The variance (ANOVA) F statistic indicated that the more overlap existed between the nitrate values in different vulnerability classes. By class difference calculation, the correct area of vulnerability assessment regions accounted for 96.4%.

The results of this study may help managers, planners, and supervisory administrations in assessing groundwater vulnerability to the pollution caused by different sources. Moreover, the results of this study help understand pollution of groundwater caused by different resources. The score scope of each parameter is a function of the region’s conditions. Hence, the score scope of each parameter in the modified DRASTICindex is determined considering land use.

The study revealed that use of the geographic information system (GIS) easily provides for preparation of DRASTICindex layers for identification of vulnerable areas. The vulnerability map is a suitable means of local and regional assessment and identification of areas prone to pollution, especially non-point pollutions for assessing the groundwater vulnerability potential in designing a supervisory network. This research also showed that the DRASTICindex could be used to prioritize vulnerable areas for identifying and preventing a higher flow of pollutants to these regions.