In this research, we employed four models DRASTIC, SI, GODS, and SINTACS to achieve the research goals and determine the level of vulnerability.
Interpretation of DRASTIC model results
In the DRASTIC method, weights have been assigned to seven parameters of groundwater depth (D), net recharge (R), aquifer environment (A), soil environment (S), topography (T), unsaturated environment (I), hydraulic conductivity (C) in GIS. Table 1 shows the weights of each of the classes of layers used in the DRASTIC method.
Table (1) ranking and weight of DRASTIC model parameters
Groundwater depth / meter
|
Range
|
0–25
|
25–39
|
39–50
|
50–63
|
63–77
|
77–92
|
92-134.5
|
Rank
|
10
|
9
|
7
|
5
|
3
|
2
|
1
|
Net Nutrition / mm
|
Range
|
3–5
|
5–7
|
7–9
|
9–11
|
11–13
|
-
|
-
|
Rank
|
1
|
3
|
5
|
8
|
10
|
-
|
-
|
Saturated zone
|
Lithology
|
Sand
|
Sandstone
|
Clay
|
Gravel
|
Shale
|
Limestone
|
Slit
|
Rank
|
8
|
8
|
1
|
10
|
3
|
3
|
1
|
Soil environment
|
Soil layer
|
Clay
|
Clay lime
|
Lime
|
Sandy lime
|
Lime silt
|
Silty
|
-
|
Rank
|
1
|
2
|
4
|
5
|
3
|
2
|
-
|
Topography %
|
Range
|
0–2
|
2–6
|
6–12
|
12–18
|
> 18
|
-
|
-
|
Rank
|
10
|
9
|
5
|
3
|
1
|
-
|
-
|
Unsaturated zone
|
Lithology
|
Sand
|
Clay
|
Gravel
|
Sandstone
|
Shale
|
Marl
|
Silt
|
Rank
|
8
|
1
|
10
|
8
|
6
|
1
|
1
|
Hydraulic conductivity, meter / day
|
Range
|
0.01–1.3
|
1.3–3.9
|
3.9–8.6
|
8.6–13
|
13-24.2
|
> 24.2
|
-
|
Rank
|
1
|
2
|
4
|
6
|
8
|
10
|
-
|
Depth of underground water (D):
The distance between the surface of the earth and the water table is called the depth of the underground water (Ahmadi et al., 2012). As the level of depth increases, the likelihood of pollution decreases (Pakhstin Rouhi et al., 2016). In this research, the groundwater depth layer was prepared using the IDW interpolation method due to its higher accuracy than other methods in ArcGIS. According to Table (1), this layer is categorized into 7 ranges from 0–25 to 92-134.5, and for each of the ranges, a proportional weight from 1 to 10 is recorded.
Net nutrition (R):
Nutrition is the amount of water penetrating the surface of the earth and entering the water table (Rahman, 2008). The intensity and passage of dissolved substances depend on the intensity and vertical movement of water into the ground (Bouwer, 1978). In this research, the Piscopo method has been used to prepare the pure feeding layer. This method considers three factors rainfall, slope (percentage), and soil permeability. To prepare the rainfall map, the average annual rainfall statistics of the stations in the region and the IDW interpolation method have been used due to greater accuracy. The slope map has been prepared using the Digital Elevation Model (DEM) and the soil permeability map has also been created based on the soil map of the study area. The ranking of maps of soil permeability, precipitation, slope, and net nutrition are presented in Table (2). According to this table, the 1st rank has been assigned to the rainfall layer, because the average rainfall of the region is 276 mm per year. Finally, based on Piscopo's equation, the prioritized maps of slope, permeability, and rainfall rank have been processed in the ArcGIS software to generate a nutrition map.
Table (2) Classification of the net nutrition layer (Piscopo, 2001)
Soil Permeability
|
Rainfall / mm
|
Slope / Percentage
|
Net Nutrition
|
Factor
|
Range / %
|
Factor
|
Range / %
|
Factor
|
Range / %
|
Factor
|
Range / %
|
5
|
High very
|
4
|
> 850
|
4
|
< 2
|
10
|
11–13
|
4
|
High
|
3
|
700–850
|
3
|
2–10
|
8
|
9–11
|
3
|
Moderate
|
2
|
500–700
|
2
|
10–33
|
5
|
7–9
|
2
|
Low
|
1
|
< 500
|
1
|
> 33
|
3
|
5–7
|
1
|
Very low
|
-
|
-
|
-
|
-
|
1
|
3–5
|
Aquifer environment (A):
This parameter is related to the characteristics of the materials configured in the saturated zone, such as porosity, type, and size of particles (Brahim, et, al. 2012). The aquifer environment and its ingredients determine the length and trend of the underground water flow system in the aquifer (Voudouris, et, al. 2010). To prepare a map of this factor, information related to the type of saturated layer has been used in 39 drilling and piping sections and geophysical and geological explorations in the study area. The aquifer environment of the study area consists of sand, clay, shale, silt, gravel, and lime and each is weighted according to the permeability which is presented in Table (1). Coarse-grained particles such as sand and gravel have a higher weight, and finer-grained particles such as clay and silt have a lower weight. In areas with higher permeability, there is no reaction between the pollutant and the soil, and the potential for pollution can be higher due to the high speed of pollutant diffusion.
Soil environment (S):
The soil layer, usually with a thickness of about 0.5 to 2 meters, has a great potential to remove and reduce the concentration of pollutants due to the relatively high microbial activity, the presence of many organic substances, and the presence of plant roots (Amir Ahmadi et al., 2019). The potential of aquifer pollution depends on soil properties such as permeability, texture and proportion of soil organic matter (Kim and Hamm, 1999). In this research, the analysis results of 39 profiles taken from the region have been used to prepare the soil map. The soil environment is defined according to its textural classification and scored based on the pollution potential. The soil environment is classified in the study range from clay with a weight of 1 to loam-sandy with a weight of 5. Given that the potential of aquifer pollution depends on soil characteristics such as permeability, texture, and percentage of soil organic matter, the presence of fine-grained materials such as silt and clay has reduced the relative permeability of the soil and limited the movement of pollutants. The coarse-grained materials such as sand have facilitated the movement of pollutants in parts of the study area.
Topography (T):
The topography factor expresses the changes in the slope of the area (Stigter. al. et, 2006). Slopes that provide a higher infiltration opportunity have a higher contamination potential (Plymale, et. al, 2002). In this research, the slope layer has been prepared by a DEM layer with a spatial resolution of 12.5 meters in ArcGIS software. The average slope between the two points is obtained by dividing the vertical distance by the horizontal distance between them.
Unsaturated zone (I):
This area is unsaturated or continuously saturated so that it can control the passage of the pollutants and their dilution (Ahmadi and Aberoumand, 2009). An unsaturated zone controls the duration and manner of the movement of pollutants and, therefore, affects the time required for the quantity and concentration of substances that come into contact with pollution (Plymale, al. et, 2002). To prepare a map of this parameter, we used the data related to the soil type of the unsaturated zone in 39 drilling and pipe-laying sections in the Mashhad plain. The method of obtaining information related to the unsaturated zone has been the same as that of the aquifer environment. However, in this case, we considered the granularity and characteristics of the sediments between the underground water level and the ground surface. Accordingly, the more permeable the unsaturated zone is, the less reaction is formed between pollutants and soil particles, and the washing speed towards the aquifer is faster and higher. Therefore, the potential for contamination has increased. According to what can be seen in Table (1), the formations of the study area are clay, marl, and silt (weight 1), shale (weight 6), sand and gravel (weight 8), and gravel formation (weight 10). So, the permeability in the areas with gravel formation is greater than in other places.
Hydraulic conductivity (C):
Hydraulic conductivity is the ability of the aquifer environment to transport water and its associated pollutants. This parameter controls the transport and dispersion of the pollutants from the injection point inside the saturated zone (Rahman, 2008). The information related to hydraulic conductivity is obtained from pumping test calculations (Khodai et al., 2015). The transfer coefficient or transferability of an aquifer layer is the amount of water that passes through a unit cross-sectional area of the aquifer layer under a certain hydraulic slope, in square meters per day (meter per day in each meter of layer thickness). In general, because in pumping tests, the value of the parameter of the water transfer capability coefficient is measured, therefore, by using the saturation thickness of the aquifer, the value of hydraulic conductivity is obtained by dividing the water transfer capability coefficient by the saturation thickness. The transmission capability map has been obtained from the results of pumping tests and using the IDW interpolation method. The saturated layer thickness map has been obtained from the interpolation of the data related to the depth of the saturated region obtained from the drilling sections using the IDW method. Hence, the hydraulic conductivity of the study area has been obtained in 6 ranges from 0.01–1.3 to more than 24.2 and with weighting from 1 to 10 according to Table (1).
Figure 4 shows each layer of the DRASTIC model and the rank of each class. After preparing each of the layers of the DRASTIC model to prepare the vulnerability map of the aquifer, the layers of the model have been combined in the ArcGIS software environment and the final map has been obtained (Figure 8). We, then, has been calculated the area of each category. Based on the results of the DRASTIC model, the studied area has been divided into five zones with very low, low, medium, high, and very high vulnerability, which are 5.81, 26.03, 44.45, 22.57, and 13.1 respectively. It includes 1% of the studied area (Table 7).
Interpretation of SI model results
In the SI (sensitivity index) method, five parameters of groundwater depth (D), net recharge (R), aquifer environment (A), topography (slope) (T) and land use (LU) have been weighted for the vulnerability of the aquifer using GIS. It should be noted that the four parameters of the saturated environment, groundwater depth, net nutrition, and topography are ranked in the SI method like the DRASTIC method, and the only difference is that in the SI method, the ranks below the criteria of these layers are multiplied by 10, that is, if the rank of one of the sub-layers in the DRASTIC method is 3, it becomes 30 in the SI method. The rank of the sub-sections of the layers of the SI method varies between zero and 100. Because its 4 parameters are the same as DRASTIC, it is avoided to describe them again and only the land use layer is described below.
Land use (LU)
The problems related to the increase of urban and industrial chemical pollutants new agricultural practices and changes in the use of land are serious threats to the environment in recent decades. Therefore, in arid and semi-arid areas where the dependence on underground water resources is greater, the destructive effect of the quality of these resources will be more intense due to the natural weakness of water and soil resources. To prepare the land use map, we have used the Landsat 8 OLI sensor satellite image, pass 159 and rows 34 and 35, related to May 14, 2021 (May 24, 1400), and with a spatial resolution of 30 meters (Fig. 1). The study area has rainfed agriculture, pastures, population centers, water agriculture, and gardens, forests, and rivers. The highest rating has been given to agricultural use because the chemical fertilizers used in the fields are washed by irrigation water and rainfall and easily penetrate the surface of the earth. The low slope in these areas also facilitates the penetration of pollutants. Residential areas, industrial facilities, and roads are in second place in terms of increasing the vulnerability potential of the aquifer. This can be a result of the pollution caused by anthropogenic factors. Also, due to its very low permeability, the river has been assigned zero rank in terms of increasing the vulnerability potential of the aquifer. The ranking of the criteria of this parameter is shown in Table (3). Table 5 shows the weight of each of the classes of layers used in the SI method.
Table (3) Land use ranking in the SI model (Sensitivity Index) (Stigter et al., 2006)
Landuse type
|
Rank
|
Agricultural areas
|
Annual crop irrigation areas (rice fields)
|
90
|
Permanent crop areas (orchards, herbal gardens)
|
70
|
Heterogeneous agricultural areas
|
50
|
Pastures and agriculture-forests
|
50
|
Anthropogenic areas
|
Industrial waste, landfill sites
|
100
|
Human-made areas (mines, shipbuilding factories, desert mines)
|
80
|
Continuous urban areas, airports, ports, highways, railways, industrial and commercial areas outside green spaces
|
75
|
Discontinuous urban areas
|
70
|
Natural areas
|
Water environments (salt marshes, salt lakes, tidal areas)
|
50
|
Natural areas (forests, semi-natural areas)
|
0
|
Aquatic areas
|
0
|
Figure (5) shows the map of each layer of the SI model and the rank of each class. After preparing each of the layers of the SI model to prepare the vulnerability map of the aquifer, the layers of the model have been combined in the ArcGIS software to create the final map (Fig. 8). Then, the area of each category has been calculated. Based on the results of the Mashhad plain aquifer sensitivity index model, created from the linear combination of SI model parameters, a vulnerability map has been prepared in five zones with very low, low, medium, high, and very high vulnerability, respectively, with 0.4, 24.63, 23.98, 18.71 and 32.25 percent of this area (Table 7).
Interpretation of GODS model results
This method uses four parameters aquifer type (G), lithological characteristics of unsaturated zone (O), depth of underground water (D), and type of surface soil texture (S). To achieve the goals of the research and prepare each of the layers used in the GODS model, we employed capabilities presented by GIS and remote sensing technology. Table (4) shows the weight of each of the classes of layers used in the GODS method.
Aquifer type (G):
This parameter evaluates the type of aquifer in terms of whether it is free or confined. The information about the type of aquifer has been obtained using the information of the water resources map. The type of aquifer is weighted according to Table 4, which includes free aquifer (weight 1), free covered (weight 0.5), semi-confined (weight 0.3), confined (weight 0.2), artesian (weight 1. 0) and without reservoir (weight 0). The type of aquifer throughout the study area is open aquifer type.
Lithological characteristics of the unsaturated zone (O):
Lithology examines the origin, structure, characteristics, and history of the solid layer of the earth studies rocks in terms of origin, age, composition, and distribution, and classifies them in terms of their physical and clinical properties. Geological maps of the studied area have been used to prepare the lithology layer. Thus, the Mashhad Plain aquifer study area includes three categories of materials in terms of lithological characteristics. the three categories are unhardened materials, hardened materials (porous rocks), and hardened materials (dense rocks). Unhardened materials (sediments) are formed by alluvial cone gravel, alluvial sand, glacial sand, windy sand, alluvial silt, and earthen soils. Hardened materials (porous rocks) are formed by chalky limestone, sandstone, siltstone, and mudstone. Hardened materials (dense rocks) are formed by karst limestone and recent lavas with igneous and metamorphic rocks. As can be seen in Table 4, each of these elements is weighted from 0.4 to 1.
The depth of the underground water level (D):
The groundwater depth layer has been produced by the IDW interpolation method due to its higher accuracy than other methods in ArcGIS software. According to Table (4), this layer is classified into 7 ranges from ˂2 to > 100 and a proportional weight from 1 to 7 has been recorded for each of the ranges.
Type of soil surface texture (S):
Using the results of the analysis of 39 profiles taken from the region, the soil environment has been determined according to its textural classification and scored based on the pollution potential. The soil environment in the study area is classified as clay (weight 1), clay loam (weight 2), silt (weight 3), silty loam (weight 5), loam (weight 7), and sandy loam (weight 9). The potential of aquifer pollution depends on soil characteristics such as permeability, texture and proportion of soil organic matter. Thus, the presence of fine-grained materials such as silt and clay can reduce the relative permeability of the soil and can also decrease the movement of pollutants. In this situation, coarse-grained materials such as sand can facilitate the movement of pollutants in parts of the study area.
Table 4
Rank and weight of GODS model parameters (PAEZ, 1990)
Aquifer type
|
No aquifer
|
Artesian
|
Confined
|
Semi-confined
|
Free and covered
|
Free
|
0
|
0.1
|
0.2
|
0.3
|
0.5
|
1
|
Lithologic characteristics of unsaturated zones
|
Intact soils
|
Silty sand
|
Aeolian sand
|
Alluvial sand
- Glacial sand and gravel
|
Coarse-grained materials (sediments) and argillaceous rocks
|
Argillaceous siltstone (sediments)
|
Limestone with gypsum
|
Siltstone
|
Sandstone
|
Limestone with gypsum
|
Indurated materials (compact rocks)
|
Igneous and metamorphic rocks
|
Present-day lavas
|
Karst limestone
|
Indurated materials (compact rocks)
|
0.4
|
0.5
|
0.6
|
0.7
|
0.8
|
0.9
|
1
|
|
Depth of groundwater level / meters
|
˃ 100
|
100 − 50
|
50 − 20
|
20 − 10
|
10 − 5
|
5 − 2
|
˂ 2
|
Texture of surface soil
|
Clay
|
Clay silt
|
Silt
|
Sandy silt
|
Sand
|
Coarse-grained gravel and sand
|
Absence of soil
|
0.5
|
0.6
|
0.8
|
0.9
|
1
|
|
Figure 6 shows the maps of each of the layers of the GODS model and the rank of each class. After preparing the GODS model layers to prepare the aquifer vulnerability map, the model layers have been combined in ArcGIS to produce the final vulnerability map (Fig. 8). Then, the area of each category has been calculated (Table 7). Based on the results of the GODS model, the study area is divided into five zones with very low vulnerability (0.93%), low vulnerability (31.11%), medium vulnerability (11.45%), high vulnerability (1.56%) and very high vulnerability (54.95%).
Interpretation of SINTACS model results
In SINTACS method is executed with seven parameters of groundwater depth (S), net nutrition (I), effect of unsaturated zone (N), soil type (T), aquifer environment (A), hydraulic conductivity (C), and topography (slope) (SV). We have used GIS and remote sensing technology to achieve the research goals and prepare each of the layers according to the SINTACS model. Given the fact that the parameters used in the SINTACS model are the same as the DRASTIC model, it is avoided to mention again how they are weighted. These two models are the same in the type of investigated parameters, but they are different in the calculation method and mathematical relationships. Table (5) shows the weight of each of the classes of layers used in the SINTACS method.
Table 5
Rank and weight of parameters in the SINTACS model (Eftekhari et al., 2019)
Depth of ground water (meter)
|
Range
|
0–3
|
3–5
|
5–7
|
7–10
|
10–13
|
13–20
|
20–30
|
30–36
|
36˂
|
Rank
|
9
|
8
|
7
|
6
|
5
|
4
|
3
|
2
|
1
|
Net nutrition (mm)
|
Range
|
0–50
|
50–100
|
100–180
|
180–250
|
250>
|
Rank
|
1
|
3
|
6
|
8
|
9
|
Unsaturated zone
|
Formation type
|
Clay, marl, silt
|
Shale
|
Sandstone and sand
|
Gravel
|
Rank
|
1
|
6
|
8
|
10
|
Soil environment
|
Soil type
|
Clay
|
Lime clay
|
Silt
|
Silty lime
|
Lime
|
Sandy lime
|
Rank
|
1
|
2
|
3
|
5
|
7
|
9
|
aquifer environment
|
Formation type
|
Clay and silt
|
Limestone
|
Sandstone and sand
|
Gravel
|
Rank
|
1
|
3
|
8
|
10
|
Hydraulic conductivity (meter/day)
|
Range
|
0.01–1.3
|
1.3–3.9
|
3.9–8.6
|
8.6–13
|
13–24
|
24>
|
Rank
|
1
|
2
|
4
|
6
|
8
|
10
|
Topography (slope %)
|
Range
|
0.2
|
2–6
|
6–12
|
12–18
|
18˂
|
Rank
|
10
|
9
|
5
|
3
|
1
|
Figure 7 shows the maps of each layer of the SINTACS model and the rank of each class. After preparing each of the layers of the SINTACS model to prepare the vulnerability map of the aquifer, the layers of the model have been combined in the ArcGIS software environment to generate the final map (Fig. 8). The area of each category is presented in Table 7. According to the SINTACS model, the Mashhad plain aquifer is divided into five vulnerability classes very low (0.44%), low (25.57%), moderate (28.58%), high (2.79%), very high (42.61%).
Validation of DRASTIC, SI, GODS, and SINTACS vulnerability maps:
To validate the vulnerability maps prepared from all four models, we have calculated the correlation coefficients between the vulnerability maps and the quality index of TDS. Information about TDS values of piezometric wells in the study area has been obtained from the Regional Water Company of Razavi Khorasan province. In this research, validation has been performed using statistical methods and calculating the correlation coefficient between vulnerability maps and TDS values in TerrSet software. In general, the concentration of TDS values of total dissolved solids is high in polluted groundwater and relatively low in drinking water (Ozler, 2003 & Mohammadi al. et, 2012). According to the classification (Todd, 1980), the TDS layer is divided into two categories of 0-1000 mg/liter class representing fresh water and 1000. 10,000 mg/liter class representing salty and polluted water. Validation results showed that all four DRASTIC, SI, GODS, and SINTACS models have high accuracy in zoning the vulnerability of Mashhad plain aquifer so that the correlation coefficients of vulnerability maps with TDS quality index are 0.996 in the DRASTIC model, and 0.995 in SI model, 0.85 in the GODS model, and 0.91 in the SINTACS model (Table 6).
Table 6
Correlation coefficient of TDS values with DRASTIC, SI, GODS, and SINTACS vulnerability maps.
Correlation coefficient
|
Model
|
0.996
|
DRASTIC
|
0.995
|
SI
|
0.85
|
GODS
|
0.91
|
SINTACS
|
Vulnerability analysis of Mashhad plain aquifer
By comparing the vulnerability classes in the four DRASTIC, SI, GODS, and SINTACS models presented in Table 7 and Fig. 9, it can be stated that the widest vulnerability in SI, GODS, and SINTACS models is related to the class of very high and in DRASTIC it is related to the class of moderate. The results of three of the four models have indicated that the class of very high vulnerability covers the largest area of the aquifer (Fig. 9).
Table 7
Area and percentage of vulnerability classes in DRASTIC, SI, GODS, and SINTACS vulnerability maps.
Model
|
DRASTIC
|
SI
|
GODS
|
SINTACS
|
Vulnerability classes
|
Area %
|
Area km2
|
Area %
|
Area km2
|
Area %
|
Area km2
|
Area %
|
Area km2
|
Very low
|
5.81
|
157.9
|
0.4
|
10.9
|
0.93
|
25.21
|
0.44
|
11.92
|
Low
|
26.03
|
707.38
|
24.63
|
669.54
|
31.11
|
840.02
|
25.57
|
690.5
|
Moderate
|
44.45
|
1207.97
|
23.98
|
651.9
|
11.45
|
309.2
|
28.58
|
771.73
|
High
|
22.57
|
613.42
|
18.71
|
508.5
|
1.56
|
42
|
2.79
|
75.38
|
Very high
|
1.13
|
30.72
|
32.25
|
876.57
|
54.95
|
1483.5
|
42.61
|
1150.4
|