Assessment of groundwater potential using multi-criteria decision analysis and geoelectrical surveying

ABSTRACT A precise map of the dispersion of the groundwater potential across each watershed can help decision-makers to exert optimal water management in each region. In this research, the potential of groundwater resources in both the Zanjanrood Catchment and the Tarom Region, located in the northwest of Iran, has been studied. Seven effective criteria including slope, land-use, drainage density, spring density, lithology, lineament density, and rainfall are considered. Criteria were first weighted using the Analytical Hierarchical Process (AHP) method and then overlaid by the Technique for Order Preferences by Similarity to Ideal (TOPSIS) model. Finally, the spatial zoning map of groundwater potential was obtained in four categories. A sensitivity analysis was performed to determine the influence of each criterion on the obtained map. The model was verified using both the spatial distribution of the high-discharged production wells and the geophysical-based geoelectric field surveys. The results indicate that the high-discharged wells (>40 l/s) in both regions are dispersed predominantly in the very good zone and, in several cases, in the good zone. Besides, the results from the two-dimensional models of resistivity and induced polarization of geoelectrical field survey are inappropriate agreement with those from the TOPSIS method. Notably, there is no suitable potential zone of groundwater in the surrounding highlands to be used in the future for drinking purposes since the highlands water supply is a strategic supply for drinking. The strategies employed in this study, the results of GIS modeling, and the geoelectrical analysis can be considered for sustainable development of similar arid and semi-arid areas since groundwater is considered as the main supplier of water in such regions.


Introduction
Recently, a water crisis has been happened in most portions of Iran with (semi) arid climates as a result of the population increase, inefficient agriculture sector, and water resources mismanagement (Madani 2014). These factors have led to severe aquifer storage depletion. If this condition continues, the water demand will not be met sufficiently in the future (Trinder and Liu 2020). To reduce the water crisis's effects on people's lives, the water resources have to be managed optimally to achieve sustainable development (Scott and Rajabifard 2017). While groundwater is the primary source of water supply in most portions of the country (Masoumi, Mousavi, and Hajeb 2021). A precise spatial map of the dispersion of the groundwater potential across each watershed could help decision-makers to conduct optimal water management in each region (Kumar, Dwivedi, and Gaur 2021). Therefore, this research tries to investigate the potential of groundwater resources in the northern portion of the Zanjan Province, Iran. The recent overexploitation of groundwater along with the resent-years droughts have led to a severe drop in groundwater levels, particularly in the Zanjan alluvial aquifer. For example, according to the reports of the Zanjan province Regional Water Company (ZRWC) in the Zanjan alluvial aquifer, a groundwater level declined as high as 34.4 m from 1991 to 2011 which is an average of 26 cm in a year. In this research, we have tried to take into account the criteria of slope, landuse, drainage density, lineament density, lithology, rainfall, and spring density to assess the groundwater potential in both the Zanjanrood Catchment and the Tarom Region. All the selected criteria except the small-discharged springs density have been widely applied for groundwater potential assessment by the previous researchers (Saraf and Choudhury 1998;Nagarajan and Singh 2009;Dar, Sankar, and Dar 2010;Magesh, Chandrasekar, and Soundranayagam 2012;Mahmoud and Alazba 2016;Sahoo et al. 2017;Sener and Davraz 2013;Arulbalaji et al. 2019;Achu, Thomas, and Reghunath 2020a). Moreover, the data availability is considered in selecting criteria. We used Analytical Hierarchical Process (AHP) and Technique for Order Preferences by Similarity to Ideal Solution (TOPSIS), respectively, to weigh and overlay the different criteria. AHP method is based on the paired comparison which quantifies the weights of decision criteria (Aju et al. 2021). It has been used in many previous researches to obtain the weight of each criterion (Sener and Davraz 2013;Rahmati et al. 2015;Razandi et al. 2015;Shekhar and Pandey 2015;Pinto et al. 2017;Achu, Reghunath, and Thomas 2020b). Besides, TOPSIS is a method based on goal programming, has been previously used to integrate data layers (Karimi et al. 2011;Rad and Busch 2011;Pazand, Hezarkhani, and Ataei 2012;Sánchez-Lozano, García-Cascales, and Lamata 2016;Zaree, Javadi, and Neshat 2019). The main advantages of TOPSIS include: (1) it considers a significant number of criteria, (2) it is easy to implement, (3) the method of determining the weights in it is optional, (4) prioritization in this method is performed based on the ideal solution, and (5) for the best and the worst criteria calculate a scalar value simultaneously (Kim, Park, and Yoon 1997). In general, data layers prepare from topography maps, geological maps, satellite images, drainage network maps, rainfall data, spring data, and Digital Elevation Model (DEM) (Ganapuram et al. 2009;Chowdhury, Jha, and Chowdary 2010;Gupta and Srivastava 2010;Preeja et al. 2011;Aju et al. 2021;Tolche 2021;Chaudhry, Kumar, and Alam 2021;Azma et al. 2021).
To assess the groundwater potential using Geospatial Information System-based (GIS-based) methods, one first must define the controlling factors which are considered as input layers to the model (Saatsaz et al. 2020). In previous studies, researchers have considered various criteria such as different effective criteria in groundwater potential, which can be due to differences in the study area and available data.
In this study, the geoelectrical field survey of resistivity and Induced Polarization (IP) in selected points has been conducted to verify and evaluate the results derived from the GIS model. The geoelectrical survey has been very successful in delineating sub-surface layering gives the resistivities, depths and thicknesses of each geoelectrical layer (Naiyeju, Oladunjoye, and Adeniran 2021). This method is considered to be the most promising and most suitable method for ground water prospecting which employed in this study.

Study area
The Zanjanrood Catchment (4857 km 2 ) and the Tarom Region (3073 km 2 ) are located in the northern portion of the Zanjan Province, north-west of Iran ( Figure 1). The average rainfall in Zanjanrood Catchment and Tarom Region is 330 and 240 mm/ yr, respectively (Meteorological Organization of Islamic Republic of Iran). The Zanjanrood River, the main drainage of the Zanjanrood Catchment, with a length of about 120 km gathers the runoff producing over the catchment and flows westwards and joins the Ghezelozan River. The Ghezelozan River with a length of about 550 km originates from Chehel-Cheshmeh spring in Kurdistan. About 165 km of this length flows through Tarom Region. Zanjanrood and Ghezelozan Rivers are two main permanent rivers in the study area (Zhao et al. 2022). Moreover, 22 sub-streams flow from the mountain in Tarom Region to the Ghezelozan River. Changing the type of rainfall from snow to rain, increasing the average temperature, and changing the temporal and spatial distribution of rainfall have caused climate change and drought in the Tarom Region so that it has reduced surface water and increased groundwater exploitation (Abdinejad 2011). Also, according to statistics from 2008 to 2010, the total annual volume of groundwater pumping in Zanjanrood catchment is 320 million m 3 . The large volume of groundwater in both regions is used for agricultural purposes (Zanjan Regional Water Organization (ZRWO)).
The rocks cropped out in the Zanjanrood Catchment are mainly Precambrian dolomitic limestone; Cenozoic upper red formations consist of shale and clay; the Quaternary alluvial deposits consist of silt, clay, and sand (Darvishzadeh 2002). The Tarom Region is a magma arc that includes volcanic and Eocene volcanic deposits, which were impressed by Oligocene intrusive masses and tolerated slow tectonics (Aghanabati (2004)). The rocks exposed in the Tarom Region at first are Cenozoic exposed consist of conglomerate and sandstone then, calcic and volcanic deposits. More volume of deposits Includes intermediate andesitic rocks (Aghanabati 2004). In the basin that is located among Tarom Mountains, there are mostly, Neogene and Quaternary alluvial deposits, unconsolidated conglomerate, and marl (Hirayama et al. 1966). Geologically, the area is located in the structural zone of Central Iran. The main geologic units of the area, in decreasing age order, are (1) Precambrian Kahar formation (Mt2), which consists silty shale and sandstone; (2) Cambrian Bayandor formation (br), which consists micaceous sandstone and shale, thin dolomite intercalations; (3) Cambrian Soltanieh formation (s1), which comprises massive dolomite with chert nodules; (4) Cambrian Barut formation (bt), which mostly consists cherty dolomite with sandstone micaceous shale and limestone; (5) Cambrian Zaigun micaceous shale formation (z); (6) Cambrian Lalun arkosic sandstone formation (l); (7) Cambrian Milla formation (Ml), which consists dolomite and dolomite limestone; (8) Permian Dorud formation (Pd), which contains quartzite, shale and sandstone; (9) Permian Routeh limestone formation (Pr); (10) Jurassic Shemshak sandstone and shale formation

Methodology
To assess a study area for groundwater potential, the primary step is to select appropriate criteria that either intensively influence both water percolation and storage or provide appropriate evidence for the presence of groundwater (Mahmoud and Alazba 2016;Oulidi et al. 2009). Therefore, considering the goal and literature review, the data and information availability, and the characteristics of the study area, we selected seven criteria, including slope, land-use, drainage density, lineament density, lithology, rainfall, and spring density. The adopted steps in the present study are shown in Figure 2. We first collected and compiled all required data to prepare the desired criteria. Then, the derived spatial layers were weighted using the AHP method. Weighted data layers were overlaid to obtain the spatial map of the groundwater potential using the TOPSIS method. After obtaining the spatial map, we performed a sensitivity analysis to highlight the role of each layer in the model results. Finally, we used both the geoelectrical field survey (Resistivity and IP) and spatial distribution of high-discharged wells to verify the model results.
The analyses were conducted using different software in each step. We have employed Expert Choice software to calculate the weights using the AHP method in this research. Then, ArcGIS10.3 and QGIS3.4 have been used to prepare data and analyze of the TOPSIS overlay model.

The criteria used to assess the groundwater potential
In this study, seven effective criteria, including slope, land-use, drainage density, lineament density, lithology, rainfall, and spring density, have been selected based on previous studies in this field and the characteristics of the study area to assess groundwater potential (Oulidi et al. 2009;Magesh, Chandrasekar, and Soundranayagam 2012;Mahmoud and Alazba 2016;Sahoo et al. 2017;Sener and Davraz 2013;Arulbalaji et al. 2019;Achu, Reghunath, and Thomas 2020b;Aju et al. 2021). In the following, the importance of these criteria, their categorization, and the rank of each category are shown in Table 1. In Table 1, rank shows the importance of classes in each criterion. The first rank, which its value is considered as 1, means the category with the highest importance in terms of groundwater potential and larger values indicate less importance, respectively. The rank of some classes is considered as half values, which means that the difference between classes is not as much as one unit in rank. Actually, in this state, the rank of the two classes is more near to each other. Consequently, fewer values specify more important classes. The reasons for the importance of each class of criteria are explained in sections 3.1.1 to 3.1.7. Furthermore, the interval and break values of the classes, in each criterion have been specified based on previous studies (Ganapuram et al. 2009;Chowdhury, Jha, and  Tolche 2021; Chaudhry, Kumar, and Alam 2021; Azma et al. 2021) and considering the characteristics of the study area. It is notable to say that the concept of rank in Table 1 differs from layer weights. The weights will be given in the next sections.

Slope
Topography slope influences the infiltration rate of rainfall into the ground and then the percolation to the groundwater (Ganapuram et al. 2009;Nistor 2019). In fact, under the same condition, in the areas with a higher slope, the infiltration rate is less. We prepared the slope layer using a DEM derived from 25 m contour lines retrieved from the 1: 25,000 topographic maps which have been prepared by the National Cartographic Center (NCC) of Iran.

Drainage density
Drainage density can also influence the infiltration rate of water into the ground, particularly during precipitation since rainfall water molecules can leave the system with a dense network of drainage quicker. In other words, the rainfall water molecules can flow with a higher velocity in the system with a dense drainage network. The development of a drainage system depends on the climate, lithology, relief, infiltration capacity, vegetation covers, and runoff intensity index. Low-permeable lithology, sparse vegetation, highrelief areas with a higher rainfall intensity lead to a denser drainage network (Bali et al. 2012). It, therefore, appears that the areas with a dense network of drainage may have a poor aquifer, particularly in the alluvial sediments. We used a 1:25,000 drainage network derived from 1:25,000 maps of the NCC of Iran to produce this layer.

Rainfall
Theoretically, the groundwater potential in the areas with a higher rainfall rate is higher than those areas with a lower rainfall rate under the same geomorphological condition. In this study, to prepare a spatial map of the rainfall, the long-run average values of annual rainfall (13 years) measured in 17 pluviometry stations scattered across the area in the different elevations were interpolated using the Kriging method (Masoumi, Rezaei, and Maleki 2019).

Lithology
The lithological character of the rocks cropped out in the area can intensively control both the infiltration and percolation of water into the underground (Mukherjee, Singh, and Mukherjee 2012). In addition, the lithology can intensively influence the storativity and permeability of the aquifers since both the porosity and fracture pattern are controlled predominantly by the lithology. The rocks/formations outcropped over the surrounding highlands of the Zanjanrood Catchment predominantly consist of volcanic, clasticvolcanic, and igneous rock, including tuff, andesite, rhyolite, granite, basalt, microdiorite, and syenite. In the Zanjanrood Catchment, the area of the plain (37.45%) which has a low slope and consists of thick alluvial and clastic sediments containing gravel, sand, silt, and clay, is larger than highlands (63.55%) with a higher slope. Note that the Zanjan aquifer which is mostly distributed around the Zanjanrood River is the common source of water for agricultural, domestic, and industrial activities in the area of interest.
In the Tarom Region, highlands (95.22%) have a higher slope where most of the outcropped rocks are Tertiary volcanic and granitoid, basalt, Porphyritic quartz, tuff, andesite rocks. The alluvial plain of the Tarom Region (4.78%) has a low slope and consists of Neogene and Quaternary deposits containing sandstone, marl, gypsum, and conglomerate. While the groundwater potential is higher in the areas with higher permeability, we tried to categorize the different formations and lithology across the area based on the table of permeability versus lithology presented by Freeze and Cherry (1979). In fact, we assigned an upper potential to the coarsegrained sediments and the karstic limestone and lower potential to the low permeable units of shale, marl, and clay sediments. The other lithological classes were categorized between them.

Lineament density
The lineaments include linear features such as fault, fractures, and geological contacts between different lithologies, which may provide a pathway for water to percolate into the ground and move through the aquifer system, particularly in the lineaments that are not filled by the soil. Therefore, the lineament density of an area can signify the groundwater potential indirectly where it supposes that the water infiltration rate increases with increasing lineament density (Magesh, Chandrasekar, and Soundranayagam 2012).
Here, to extract the lineaments across the area we used band-8 of Landsat8 satellite images which is a panchromatic band with 15 m spatial resolution. Directional filters are used to edge detection in satellite images. To obtain the directional filter kernels we used the forward (Equation 1), backward (Equation 2), and central (Equation 3) numerical differentiation as follows.
( 3) where f 0 x ð Þ is the function value in a pixel, f x þ 1 ð Þ presents the function value in the next pixel, f x À 1 ð Þ is the function value in the previous pixel, and Δx signifies the distance difference between two identical points in each pixel.
In directional filters, the pixel values change depending on the weight of elements in such a manner that the features situate in filtered directions become more visible. Therefore, at the first, the image filtered at angles of 0, 45, 90, and 315 degrees. Afterward, lineaments were obtained by automatic lineaments extraction. The results were visually corrected and adapted to the specifications of the area. Finally, the lineament density map was obtained. We also verified the lineament network obtained from the satellite with lithology. Using the map of the lineaments overlaid on the geological map, it is observed that most of the derived lineaments are scattered on the hard-rock outcrops, not the soft rocks (i.e. alluvial sediments and marl) implying that the model is working correctly.

Land-use
The type of land-use is effective in water penetration into the ground (doer et al. 2006;Masoumi, Coello Coello, and Mansourian 2020). For example, due to plowing in the agricultural land, the water infiltration rate is higher than residential areas covered by either asphalt or pavement (Osuji et al. 2010;Pingping et al. 2013). The land-use type in the study area in accordance with the order of infiltration includes Irrigated agriculture, Mixed cropland and orchards, Dryland farming, Mixed dryland farming and pasture, Low dense pastures, Medium dense pastures, High dense pastures, Mixed forest land, and Mixed urban or builtup land. These categories were extracted based on the United States Geological Survey (USGS) (Anderson 1976) land-use/land cover classification system.
The land-use layer was obtained from the Management and Planning Organization (MPO) of Zanjan Province. The boundaries in the layer were then visually evaluated using Landsat8 satellite images. Next, the land-use classes changed to the USGS standard classifications. Finally, the layer was entered into the GIS model after the necessary edits.

Spring density
The last used criterion is the density of the smalldischarged springs across the area. Over most portions of the highlands in the area, a large number of smalldischarged springs (less than 5 l/s) exist. Totally there are 6028 and 977 small-discharged springs in the Zanjan Catchment and Tarom Region, respectively. There is only one spring with the discharge of 15 l/s in the area while the rest springs are all less than 5 l/s. This issue allows us to take the spring density into account as an effective layer in assessing the groundwater potential. In fact, when the system has a low vertical permeability, the water cannot percolate to the deeper aquifer layers and then a good aquifer does not form. Therefore, water infiltrated into the shallow depth prefer to discharge as small springs around the valleys and then cannot form a large integrated aquifer throughout an area. Therefore, it could be concluded that a rich aquifer may not be formed in those areas where there are a larger number of small-discharged springs. Note that the presence of large-discharged springs such as happen in karst terrains is evidence of the presence of a rich aquifer with high permeability. Note that in our area the karstic lithological layer (i.e. limestone and dolomite) has no wide outcrop. Here, to calculate the spring density, by considering a neighborhood with a radius of 3500 m around the raster cell, the number of springs located in this neighborhood is divided by its area.

Weighting the criteria
In the present paper, the AHP method was used to assign the weight to each data layer. AHP is one of the most popular methods to determine criteria weights in the Multi-Criteria Decision Making (MCDM) (Ujoh, Igbawua, and Paul 2019;Masoumi 2022). The first step in the AHP method is to create the hierarchical structure of the reviewed subject to indicate the goal, criteria, and the relationship between them (Vargas 1990;Goepel 2018;Sevinc, Gür, and Eren 2018). In the AHP process, at first, the pairwise comparison of the alternatives was conducted using the designed questionnaires according to AHP's common questionnaires. Judgments of various specialists related to the groundwater potential such as hydrologists, hydrogeologists, environmentalists, and GIS were gathered in this research. Then, the pairwise comparison matrix (A) was created based on Equation (4), where a ij in the A matrix is the preference of criterion i over criterion j.
In Equation (4), A is the pairwise comparison matrix and a ij s show the elements of the comparison matrix, and n indicates the number of criteria. In AHP, experts usually specify the importance of criterion i over j (a ij s) as values of 1, 3, 5, 7, 9 to represent the degree/intensity of importance and the values of 2, 4, 6, 8 as median values between two degrees. (Saatsaz et al. 2020). Next, the vector of weights,w = [w 1 ,w 2 ,w 3 , . . ., w n ] were calculated based on Saaty's eigenvector method. Finally, the consistency of judgments of different experts were examined and finalized. The Consistency Ratio (CR) shows how consistent the judgments are. If the CR is much in excess of 0.1 the judgments are untrustworthy because they are too close for comfort to randomness and the exercise must be repeated. To calculate the CR firstly the λ max (the highest eigenvalue of the matrix) has to be calculated as Equation (5) (Saaty 1980).
where m represents the number of independent rows of the matrix, A represents pair-wise comparison matrix, and v means the matrix eigenvector.
Then, the Consistency Index (CI) can be calculated as Equation (6).
If the matrix is perfectly consistent then CI = 0. When dealing with rising number of pair-wise comparisons the possibility of consistency error is also increasing. Thus Saaty (1980) suggested another measure the CR (Consistency Ratio) that can be calculated as Equation (7).
where RI is the random index whose values for m number of evaluation criteria are given in Table 2 according to Saaty (1980). The suggested value of the CR should be no higher than 0.1 (Franek and Kresta 2014).

Overlay analysis
To overlay the weighted data layers, we used the TOPSIS model which is widely used as a multi-criteria decision-making method to integrate data layers (Kumar and Garg 2018). In this method, the multicriteria discontinuous analysis is initially assessed the m alternatives by the n criteria and then the alternatives are ranked in accordance to the similarity to the ideal solution (these ranks are identified in Table 1). This technique is based on the fact that the selected alternative has to be at the least distance from the ideal solution (the best possible condition) and the farthest distance from the negative ideal solution (the worst possible condition) (Lin 2010). To overlay the layers, they first have to be converted to the raster format. For detailed information about the methodology, one can refer to Boran et al. (2009), Abedi and Norouzi (2016), and Sánchez-Lozano, García-Cascales, and Lamata (2016). Then, the following steps will have performed: (1) The decision matrix including m alternatives and n criteria are created.
r ij ¼ x ij ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi P m i¼1 x ij q ; i ¼ 1; 2; . . . ; n; j ¼ 1; 2; . . . ; n (8) where r ij is normalized value for each criterion and x ij signifies the numeric value derived from i th pixel based on j th criterion. (3) The weighted normalized decision matrix is calculated using derived weights from AHP. (4) Both the positive and negative ideal solutions are estimated using Equation (9) (5) To determine the distance of each pixel from the positive and negative ideas, Equation (10) where v ij is arrays of the weighted matrix, v þ i signifies the alternative of the positive ideal, v À i is the alternative of the negative ideal, and d þ i , d À i are the distance from positive and negative ideals, respectively. Note that the positive and negative ideals are the best and worst places for each sub-criterion, respectively.
The relative distance (C) is calculated using Equation (11).
where C is the proximity value of the i th pixel to the ideal solution. The value of C changes between 0 and 1. Finally, prioritizing of the criteria is performed based on the C-level where C = 1 represents the highest rank and C = 0 signifies the least rank.

Sensitivity analysis
Sensitivity Analysis (SA) tries to assess the dependence of model output to the weights of each input parameter (Chen, Yu, and Khan 2010;Saatsaz et al. 2020).
In the SA, the weight of each criterion was regularly changed by considering a total sum equal to 1. Then, for each weight change, the model has been recalculated to obtain the effect of each layer on modeling results. When the final model is obtained, the SA is conducted by changing the weights of the parameters and exploring their effects on the answer. The one-at-a-time (OAT) method is a common approach for SA, which changes weights of each factor one by one to investigate the effect on the outputs (Masoumi 2022). All other criteria should be considered constant for the comparability of the results. Changing the weight of criterion i is defined as Equation (12).
where W i PCL denotes the alerted weight of the i th criterion based on the Percent Change Level (PCL) of weight and W i shows the main weight achieved by the AHP method for the same criterion. PCL defines the variation in the weights which are considered as ±10 and ±20 in this research.
Since the sum of all criteria weights at any PCL should always be equal to 1, when weight is changed, the other criteria are adjusted proportionally to satisfy this constraint. So, the weights of the other criteria can be computed as Equation (13): where W j PCL is the new weight of the j th criterion after changing the i th criterion based on PCL (Chen, Yu, and Khan 2010).
Lastly, those parameters that have the greatest impact on the final model are recognized to be carefully re-prepared. Finally, the last answer is obtained using the re-prepared criteria (data layers).

Results and discussion
The criteria weights have been obtained by performing the AHP method for both areas. The weights of criteria and consistency ratio obtained for both areas are given in Table 3. A consistency ratio of less than 0.1 signifies the absence of disturbance in the pair comparison calculations.
As can be seen in Table 3, the lithology has a higher weight in the Tarom Region than the Zanjanrood Catchment because according to the extended area, the Tarom Region has a higher lithological diversity. Due to the roughness of the Tarom Region and the steep slope, the slope has a higher weight in this region. The land-use in the Zanjanrood Catchment is more diverse and this can be the reason for higher weight. The lineament density in both areas has the lowest weight due to the lack of karst development and extensive soil cover.
The effect of weight on data layers is shown in Figure 3. Then, the weighted layers were overlaid using the TOPSIS method to obtain the spatial groundwater potential map. The groundwater potential maps derived by the TOPSIS method for both Zanjanrood Catchment and Tarom Region are shown in Figure 4.
The groundwater potential maps for the Zanjanrood Catchment and the Tarom Region were classified in four categories including very good, good, weak, and very weak (Figure 4). In the map of the Zanjanrood Catchment, the area of potential zones of very good, good, weak, and very weak are equal to 215 (43.4%), 1549 (31.91%), 1678 (34.57%), and 1412 km 2 (29%), respectively. The zone of very good potential covers the small portion of the Zanjanrood Catchment. This zone is mostly dispersed around the Zanjanrood River were consists of coarse-grained sediments and is dominantly recharged by the river. Figure 5 shows the lithological map of two regions for more information. The zone of good potential covers a large area of the Zanjan Plain somewhere between the Zanjanrood River and the hillside of the surrounding highlands, particularly in the southeastern and middle portions of the plain. These portions seem to receive the mountain front recharge coming from the surrounding highlands. Notably, the highlands of the Zanjanrood Catchment are not categorized as areas with good or very good groundwater potential. This may be due to that (1) the solid is filled the fracture and joints since there is a relatively thick soil covers most portions of the highlands, (2) the highlands has a steep slope, (3) the absence of widespread outcrop of the karstic carbonates rocks, and (4) lack of long extension of the fractures trough the deeper parts.
In the groundwater potential map obtained for the Tarom Region, the zone with a very good potential covers a very small portion (0.26%) of the total area, particularly around the Ghezelozan river. This also most probably is due to the continuous recharge from the river and the presence of the coarse-grained sediment which are settlement by the river during a recent geological time scale. The zone with a very poor potential covers most portion of the region (91.71%) implying that the groundwater potential is generally weak in the area. This zone corresponds with high-slope nonkarstic highlands discharging numerous smalldischarged springs suggesting a media with low permeability. Totally, the zones with good and weak potentials cover a small area. The very weak zone covers a wide range of the catchment, which is consistent with highslope non-karstic highlands and numerous small springs.
SA was performed to evaluate the sensitivity of the TOPSIS model to changes in the weights of subcriteria. The results obtained from the SA in the Zanjanrood Catchment are presented in Figure 6. The results indicate that the criteria of lithology and small-discharged spring density have the highest effect on the model while the slope has the lowest influence on the model. However, the model sensitivity to the changes in weights of criteria decreases from the lithology to slope criterion as lithology (0/7497) >spring density (−0/6516) >drainage density (−0/ 0152) > rainfall (0/0152) > lineament density (0/ 0047) > land-use (0/0004) > slope (−0/0003).
.This different pattern in both regions indicates that the interaction of the different layers on the model is complicated under the same modeling procedure. This may be due to that, the most portion of the area in the Zanjanrood Catchment is covered by the alluvial deposits with a coarse-grained size while in the Tarom Region, the alluvial sediments with a coarsegrained size are less. The area of the hard-rock outcrops which have a denser network of fractures and a steeper slope in the Tarom Region is greater than that in the Zanjanrood Catchment.
To verify the groundwater potential maps derived from the TOPSIS, we used both the spatial distribution of the large-volume production wells and the geoelectric methods (resistivity and IP). While the presences of a large-volume production well demonstrate the presence of a rich groundwater potential, we overlaid the spatial distribution of the large-volume production wells (Q > 40 l/s) on the derived maps of the groundwater potential in both study areas, all the large-volume wells are interestingly dispersed within both the zones with very good and good potentials in the Zanjanrood Catchment and Tarom Region as can be seen in Figure  8. Noteworthy, none of the wells are dispersed in the very weak potential zone. To test the results of the TOPSIS model, we also used the geoelectric methods of resistivity and IP. These methods are commonly used in the investigation of near-surface resources such as groundwater and mineral exploration via transferring the electric current into the ground with an artificial source and then measuring the potential differences (Kearey, Brooks, and Hill 2013;Srigutomo and Trimadona 2016;Attwa and Ali 2018;Hasan et al. 2020). Here, we conducted a geoelectrical field survey of resistivity and IP along a profile perpendicular to the boundary of zones with very good and very poor groundwater potentials obtained from the TOPSIS model for the Zanjanrood Catchment as shown in Figure 9.
Surveying of resistivity and induced polarization were performed by Dipole-Dipole and Pole-Dipole geoelectrical arrays along with a profile with a length of 460 m and 10 m electrode distance. Twodimensional inversion modeling of resistivity and induced polarization were implemented using software RES2DINV. A combination of both the dipoledipole and pole-dipole arrays was used to increase both the lateral and vertical resolution. The derived vertical cross-section from the two-dimensional inversion models of resistivity and IP are shown in Figure 10.
Theoretically, an aquifer with fresh water (i.e. Electrical Conductivity (EC) of <1000 µS/cm) usually shows a resistivity ranging from about 32 to about 256 Ωm (Adepelumi et al. 2009). If the EC of water increases, the resistivity will decrease since for saline water (EC>10,000 µS/cm), it may decrease to less than 10 Ωm (Breier, Breier, and Edmonds 2005) which this range coincides with the resistivity of marl and clay units. To separate the saline water aquifers from the marl and clay layers, IP data is commonly used since usually high IP values are related to marl and clay units, not saline aquifers. However, we considered a resistivity ranging from 10 to about 30 Ωm as an aquifer in the area which is in accordance with the EC value of about 5600 µS/cm related to the pumping water from the production wells nearby the geoelectric profile. Therefore, we took into account the region corresponding with 10-30 Ωm resistivity values as an aquifer system, a region signifies with black boundary in Figure 10. Notably, the aquifer is completely located within the zone with a very good groundwater potential obtained from the TOPSIS model implying that the TOPSIS model is reasonable. The aquifer system is dispersed in a shallow depth ranging from about 5 down to 20 m near the Zanjanrood river. The depth of the water table corresponds with the measured depth in the production wells (about 6 m). There is also a good agreement between the coarse-grained sediments dispersed around the river and the location of the aquifer. Moreover, the production wells used for agricultural purposes over the area are only dispersed over an area where is distinguished as an aquifer region by the geoelectrical survey. In a zone with very good groundwater potential (Figure 10), the resistivity values are less than 10 Ωm where the chargeability values increase up to about 25, particularly in deeper parts. We believed that this zone consists of marl and clay units in which the aquifer cannot be formed since (1) in accordance with the field observations, this zone consists of fine-grained units (Marl), and (2) there is no production well drilled.

Conclusions
This research attempts to assess the groundwater potential using RS data and GIS spatial analysis in both the Zanjanrood Catchment and the Tarom Region, the northern portion of the Zanjan Province in the Northwest of Iran. The groundwater potential has been investigated using seven effective criteria, including slope, land-use, drainage density, lineament density, lithology, rainfall, and spring density. The criteria were first weighed using the AHP method and then overlaid by the TOPSIS model. Derived groundwater potential maps for both the regions were classified into four categories of very good, good, poor, and very poor. Overall, the results suggest the lack of rich aquifer      In addition, the zone of very good groundwater potential covers the smallest portion of both the regions, particularly near rivers that flow through the alluvial plains. The good zone in the Zanjanrood Catchment is located between the highlands and plain. This zone covers a large portion and plays an important role in the area. However, the good zone in the Tarom Region covers a small portion. The poor zone covers the largest portion of the Zanjanrood Catchment due to fine grain geology layers that are not a suitable medium to form an aquifer system. Reasons including steep slope in highlands, lack of extent karst formation, relatively thick soil cover are caused expansion of very poor zone in both of area. Derived models from TOPSIS were also tested via using both the spatial distribution of the large-volume production wells scattered over the alluvial sediments and the geoelectrical survey of resistivity and induced polarization using Dipole-Dipole and Pole-Dipole arrays. Large-volume production wells in the Zanjanrood Catchment are approximately drilled within the zones with good and very good groundwater potential categories. 29.62% of wells are located in a very good zone, 54.32% in a good zone, 14.81% in a poor zone, and 1.23% in a very poor zone. In the Tarom Region, there are only three large-volume wells which all are dispersed within the zone with a good groundwater potential. Furthermore, the aquifer zone identified by the geoelectrical survey is in complete agreement with the zone with a very good groundwater potential derived from the TOPSIS model.
In general, this study shows that the groundwater potential in both the regions is not high enough since the zone with a very good potential covers a very small portion of the regions corresponds with the current alluvial aquifers, which are under over extracting and there is no rich aquifer developed in the highlands. Therefore, optimal management of the current aquifers and similar regions is necessary to meet the water demand in the area of interest for the future and to prevent the further depletion of aquifers and sustainable development. The methodology of GIS modeling in this study and the geoelectrical analysis can be considered for sustainable development of similar arid and semi-arid regions since groundwater is considered as the main source of water demands in such areas.

Data Availability Statement
Due to privacy and ethical concerns, neither the data nor the source of the data can be made available.

Disclosure statement
No potential conflict of interest was reported by the authors.

Notes on contributors
Marzieh Shabani received master's degree in Geophysics from Institute for Advanced Studies in Basic Sciences, Zanjan, Iran in 2017. She is currently working toward the PhD degree.

Zohreh Masoumi received her PhD degree in Geospatial
Information Science (GIS) from Department of Geomatic engineering, Khajeh Nasir University of Technology (KNTU), Tehran, Iran. Her research to date has been on the fields of disaster management, resource management, artificial intelligence and its applications, and the production and processing of spatial data. She is currently an Assistant Professor in department of Earth Sciences in Institute for Advanced Studies in Basic Sciences (IASBS), Zanjan, Iran. She has authored over 25 peer-reviewed articles in the international journals. She is also an editorial board member and reviewer of some international journals in the fields of geo-spatial information science and remote sensing.
Abolfazl Rezaei received a PhD degree in Hydrogeology from Shiraz University, Shiraz, Iran in 2014. He is currently an Assistant Professor in the Department of Earth Sciences, Institute for Advanced Studies in Basic Sciences (IASBS), Zanjan, Iran. He has authored around 20 peer-reviewed articles in international journals. His research interests include the study of physical characteristics of groundwater systems, large-scale climate variability impacts on water resources, and hydrogeology of karstic aquifers.