Introduction

Growth of population and unbridled expansion of cities, together with global warming, reducing running water resources, boundless biological pollution, diminishing phreatic zones, have drawn attentions of experts and managements to karstic formations and formations waters. These karstic formations cover 11 percent of Iran, a point of much importance in studying water resources. Studying karstic resources follows the management and maintenance goals of these resources. Therefore, identifying these resources is an ineluctable necessity for purposes of consumption management (Aguilera et al. 2009; Allocca et al. 2014). Considering the importance of volume and process of recharge in development of karst, evaluation of recharge was the primary way for identifying potential locations (Allocca et al. 2014).

The recharge of karst aquifers in terms of quantity and spatial distribution depends on various factors such as climate, vegetation. The choice of a suitable method for assessing the amount of water infiltrated into karst areas is debated among researchers. Multi-parameter methods are being developed using the geographic information system tool. Groundwater recharge is an issue that has been systematically raised in various journals, especially since the mid-1980s (Healy 2010). The recharge of karst aquifers in southern Spain has been evaluated by the APLIS method (Andreo et al. 2008). In 2009, groundwater potential in northern Jordan was based on layers of altitude, soil, fault density, geomorphology, and geology using the AHP method and divided the area into 5 classes “very bad, bad, average, good, and very good” (Awawdeh et al. 2014). In the same year, aerial photographs, geological maps, geomorphology, soil hydrogeological groups, land use and vegetation, and drainage network maps were used to determine the potential of groundwater resources in a region in India and the area was divided into three areas: good, medium, and poor (Nagarajan and Singh 2009). In 2010, a multi-criteria decision-making and GIS method was used to determine the appropriate groundwater supply areas in Tunisia (Chenini et al. 2010). In Sri Lanka, the method of geographic information system was used to determine the potential of groundwater nutrition (Senanayake et al. 2016). In the same year, studies were conducted to determine the potential of groundwater in a region of India using multi-criteria decision-making MCDM and AHP method (Kumar et al. 2016). Studies show that methods based on geographic information systems are suitable for producing groundwater potential maps, especially in areas where we lack data (Russo et al. 2015). The results of all these studies have been published in reputable scientific journals.

In geological formations, a usual measuring of recharge is possible by time series in which level of underground water is measured (Bouraoui et al. 1998); heterogeneous karstic rocks, however, make this method unreliable in karstic regions, and therefore, potential locations for development of karst formations cannot be specified (Burke 1995). In measuring recharge rate, conventional methods (evaporation and transpiration, natural isotopic trackers, chemical trackers, artificial trackers, dividing precipitation by time) or numerical methods encounter limitations and problems such as lack of accurate periodic data (Conrad et al. 2004). The key to this problem lies in using methods based on GIS. This method uses variables such as geological characteristics, altitude, slope, lithology, and fracture patterns and is often employed in explaining spatio-temporal distribution of karst recharge (Bakalowicz 2005; Carrasco et al. 2014). Thanks to this method, many advances have been resized in specifying potential regions for karstic development and many researchers have used it for measuring recharge (Duran et al. 2004; Farfán et al. 2010; Fetter 2000; Galvão et al. 2015; Gerner et al. 2012; Hartmann et al. 2014; Jirkama 2007; Kirn et al. 2017; Lee et al. 2006; Marechal et al. 2006). Offering a more accurate and faster analysis, this method facilitates measuring recharge rate of underground water and development of karstification (Kirn et al. 2017). More particularly, IGME (Instituto Geologicoly Minero de Espana), a Spanish geological institute, has developed APLIS, a new analysis for measuring real recharge rate in karstic formations and accurate potential development patterns by considering natural characteristics of the aquifer (altitude, slope, lithology, infiltration, soil) (Carrasco 2005). Unlike the established methods, this method is not dependent on accurate, periodic data, a fact which makes it much economical in terms of water resource management (Martos-Rosillo et al. 2008). In addition, this method allows for mapping the spatial distribution of natural recharge in the aquifers and development of karstic formations (Martos-Rosillo et al. 2013). APLIS also can perform accurate examination of natural recharge in carbonate minerals (including direct recharge, local concentrated recharge via shallow water holes, and even indirect recharge from the substrate surface waters) (Mejías et al. 2012). APLIS can be regarded as a novel method among numerical formula for calculating recharge that can be called calculation method for aquifer recharge by remote sensing and GIS (Díaz-Guanche et al. 2013). The problems associated with conventional methods are no longer the case here thanks to assessment that considers natural properties of a given karstic aquifer (altitude, slope, lithology, gradual infiltration, fracture, and soil) to evaluate infiltration. APLIS was used in 8 aquifers in Spain. In this trial experiment, results matched actual recharge rates (Moussavi-Harami et al. 1992). The method was tested in a various place worldwide with success: SE Spain (Pedrera et al. 2017 Ramazani Oomali et al. 2008), Greece (Saadat 2014), Cuba (Samper et al. 2005), NE Oman (Scanlon et al. 2002), Mediterranean (Andreo et al. 2008), and southern Spain (Scibek and Allen 2006; Andreo et al. 2015).

This research also examines a spatial identification of regions for developing karsts in Mashhad on a 1:250,000 scale map. The local specifications and lack of periodic data (due to recent severe reduction of water resources) make APLIS a necessary method. The importance of this study is particularly tangible for provision of water resources for drinking and agriculture in the region.

Area under study

General features

The area is situated between eastern latitudes of 58°30ʹ and 60° and northern longitudes of 36° and 37°. It lies in drainage basin of Kashfrud which flows through Binalud and Hezar Masjed mountain ranges (Fig. 1). In topological terms, the area includes both plains and heights, with lowest spot standing 462 m AMSL. Heights are vast, including two distinct ranges with highest point standing 3310 m AMSL. City of Mashhad lies at center of the area under study and is located at 1050 m AMSL. Climatologically, the area is affected by different air masses, giving it a unique climate. Mashhad-Neyshabur plain displays cold, arid climate, while Mashhad-Quchan plain is characterized by cold, semi-arid climate or steppe climate. A small portion of Binalud and Hezar Masjed mountain ranges represents cold, wet climate, while the whole area under study is characterized by changing but mild climate, tending toward cold and arid conditions and by dry, hot summers and cold, wet winters. Direction of winds is mostly from SE to NW. Maximum temperature during summer is 43 °C, while winters can be as cold as  − 23° C.

Fig. 1
figure 1

Geographical coordinates of the area under study

Full size image

Geological features and structure

Geologically, the area is characterized by deposits from Alborz zone and central Iran (west to the area under study) in Paleozoic period. From the advent of Jurassic sea to the end of Cretaceous period, deposition took place in Kopedagh (Shabanian et al. 2010) (Fig. 2). Stratigraphically, there existed a shallow lake in this region in Iran from The Liassic Epoch till early Oligocene Epoch.

Fig. 2
figure 2

Geological map of Mashhad: 1:250,000

Full size image

Thick consequent marine depositions in the basin extended to 6000 m and Jurassic sediments are conformable to the Oligocene are facing each other. The marine deposits extend to Upper Cretaceous, and only mild Epeirogenesis movements occurred toward the end of the Jurassic and beginning of the Cretaceous. No significant gap has been found in depositing processes in this region, particularly in deeper deposits. The Upper Jurassic units (Mazduran formations) are the main origin of gaseous hydrocarbons. With retreating seas in the Paleocene, much of the region emerged but no considerable folded structures are detected. A regional land subsidence accounts for the reflow of sea into Kopedagh. Structurally, compressive stress stretches from N-NE to S-SE which primarily causes shortening and folding. With later increase in pressure, trust faults have been created in some place (Shabanian et al. 2010). Most of these faults are transcurrent and are characterized by reverse component, and sometimes, normal structures although some thrust faults are also found to have folded structures and are considerably older than transcurrent faults.

Materials and methods

Many methods are available for identifying karstic regions and for calculating recharge rate. The local specifications and lack of periodic data (due to recent severe reduction of water resources) make APLIS a necessary method. However, a systematic identification of determinants in form of base maps, one that aims to maintain and manage consumption of these resources, will prove less time-consuming and more economical. The key to such identification is drawing base patterns in form of base maps which would use APLIS to consider geological, lithological, pedological, tectonic, and topological properties of the area. The mathematical formula in this method involves a combinatory or adaptive layer that displays main infiltration and development parameters of karst in a certain area, because karstic development depends on a host of physical properties of the aquifer although these variables vary in terms of weight and importance in recharge. This means that information layer for each variable is assigned a score according to APLIS ranking system and are then combined using APLIS formula. The different multi-principle formulas (such as regression analysis (linear correlation with least squares), ideal point analysis, weighted linear combination) are used to determine the weight of each variable. Recharge rate and, consequently, the development of karst in formations are calculated by combining natural specification of each aquifer: AMSL altitude, slope, lithology, gradual infiltration, and soil. After inserting information layer of each variable into GIS, mean recharge rate and development of karst are determined by APLIS formula (Bakalowicz 2005):

$$R = (A + P + (3 \times L) + (2 \times 1) + S)/0.9$$
(1)

where R is mean recharge rate, A is altitude, P is slope, L is lithology, I is preferential infiltration, and S defines the soil type. The weight of each variable represents its importance in determining recharge rate (Widiastuti 2012). Thus, lithology is three times more important than altitude (AMSL), slope, and soil, while gradual infiltration is two times more important than altitude, slope, and soil.

Results and discussion

Conversion of an issue into a hierarchical structure stands as the most important part of analytical hierarchy process. In identifying karstic regions, our goal is to specify water resources by considering infiltration percentage or recharge in formation. Criteria and sub-criteria include those variables which make options differ. Each variable is characterized by a certain weight, and every region is credited according to criteria. In this way, mean recharge rate (R) for each aquifer (in percent) is used to identify karstic regions by considering mean R for each spatial unit on the map. There are ten ranks for each variable, and each is assigned a number between 1 and 10 in form of an arithmetic progression (Bakalowicz 2005). 1 and 10 indicate minimum and maximum infiltration, respectively. Therefore, in order to determine the penetration (R) in the APLIS equation, firstly, the information layers in the APLIS equation should include: elevation (A), slope (P), lithology (L), preferential penetration (I), and soil type (S) in the GIS environment. Then, the appropriate valuation for each layer of information must be done for use in the APLIS equation. The basis of the valuation is based on the factors influencing the penetration. The following factors have been investigated:

Altitude (A)

This variable is divided into units of 300 MS (Table 1), and each unit is indicated by a score between 1 and 10 (Bakalowicz 2005). Minimum altitude is 462 m; therefore, this variable extends through 6 units (5–10). Given the high altitude of the bulk of the area under study, it is generally explained by large numbers in the rank, with 10 being the frequent indicator (Fig. 3). The layer for this variable and scores for the units show that higher altitudes receive more precipitation, meaning that recharge happens in a larger scale and more karstic formations are achievable. Also, it is understood that no significant difference is visible among units higher than 3100 m AMSL.

Table 1 The altitude variable and scoring
Full size table
Fig. 3
figure 3

Altitude map and scoring elevation layers

Full size image

Slope (P)

Using GIS software, slope map was produced by inserting DEM map of altitude. The result was ranked in 9 units (Table 2). Score 7 was dropped in scoring this variable. Ranking this variable in varying distance blocks follows the ranking used in Andles Biological Data System where slope is divided into 9 ranks and score 7 is omitted (Bakalowicz 2005). Assigned scores become smaller as slope decreases. In other words, sharper slopes are associated with smaller recharge rate (Fig. 4). Therefore, minimum slope is represented with largest score (equal 10). In high altitudes we see 3, 4, and 5, scarcely getting any larger.

Table 2 The slope variable and scoring
Full size table
Fig. 4
figure 4

Slope map and scoring

Full size image

Lithology (L)

Karst development often occurs in calcareous and dolomite carbonate formations as well as in evaporation rocks. The existence of any factor such as joint, slot, and fault that causes more water to penetrate these hard formations will increase the development of karst in each formation. But in the APLIS method, any rock that can penetrate is capable of karsting, which is one of the drawbacks of this method. Field observations as well as stratigraphical studies let us use lithological and tectonic properties (Table 3) of the region under study on a scoring scale of 1 to 5 (Fig. 5). Very few fractures in the area are responsible for reducing the effectiveness of this variable,

Table 3 The lithology variable and scoring of this variable
Full size table
Fig. 5
figure 5

Lithology map and values of the layers

Full size image

Soil (S)

Table 4 shows kinds of soil and the scoring system. Evolution of soil means thicker formations and smaller grains, hence, less infiltration. Thus and according to field observations, scoring matched the soil typology (Fig. 6). In fact, the high soil thickness due to its evolution on the surface of the earth and the fine grain components of the soil cause delays in the penetration process and consequently evaporation increases.

Table 4 The soil variable and scoring of this variable
Full size table
Fig. 6
figure 6

Soil map and soil values

Full size image

Preferential infiltration (I)

This layer is the sum of scores for three layers: slope (S), fractures (F), and cover (C).

$$I = C + F + S$$
(2)

The presence of any joint, fracture, and fault increases the penetration of hardened formations and increases the karst formation of the desired formation. These three variables are scored according to their contribution to gradual infiltration (Table 5). Fractures in three ranks are scored according to the distance from the fracture. This was done with a buffer tool in the GIS environment. Slope was also scored in 9 units where each means a percentage fraction of the slope. Increasing slope increases the speed of runoff and decreases penetration time and decreases penetration. In case of cover, rocky areas got 1, while soil areas got 0. In fact, in determining the factor of penetration preferential, the soil factor is an inhibitory factor in penetration and increasing time for evaporation of water. Eventually, the mean value of these three layers was used to determine the gradual infiltration of various areas (Fig. 7). In this way, it was found that gradual infiltration occurs in regions characterized by limestone and mild slopes, making them prone to form karstic structures. In this scoring system, 1 and 10 represent minimum and maximum gradual infiltration, respectively.

Table 5 The preferential infiltration variable and scoring of this variable
Full size table
Fig. 7
figure 7

Preferential infiltration map and values

Full size image

Calculating recharge rate (R)

In order to measure recharge rate (percentage) of existing aquifers and, then, to identify potential karstic developments, the information layers were integrated by APLIS method in GIS (Fig. 8). For this purpose, the information elevations (A), slope (P), lithology (L), preferential penetration (I), and soil type (S) are evaluated as a raster, according to the APLIS equation (R = (A + P + (3 × L) + (2 × I) + S) /0.9) combined in the RASTER calculator tool, and the penetration rate (R) is obtained in percentage terms.

Fig. 8
figure 8

Approximate recharge rate in APLIS in the area under study

Full size image

According to our calculations, minimum, maximum, and mean recharge rate in our study are 15, 96, and 49%, respectively. In terms of location, the mean recharge rate in our study is 49%. This means that it falls between 32% and 54%, the acceptable range suggested by other researches (Bakalowicz 2005).

Conclusion

In order to locate karsts, the chief determinants were studied in form of digital data in GIS on a map of locations. The results determine the recharge percentage and karstic regions on a 1:250,000 map of Mashhad in NE Province of North Khorasan. The applied method was APLIS hierarchical analysis. Accordingly, parameters such as altitude, slope, lithology, soil, and gradual infiltration (distance from fracture, slope, and presence/absence of soil) were assumed to be influential factors in recharge and development of karsts. The results were classified into recharge (5 classes) and karstic development (5 classes). Considering the correlation between development of karstification and infiltration into formation, we can see that development of karstification means increase in infiltration, while reduction of karstification means less infiltration. Therefore, mean recharge rate (R) in percentage for each spatial unit helps us identify locations with potential for becoming karstification (Willaarts et al.2012) (Table 6).

Table 6 Qualitative classification of recharge and spatial development of karstification
Full size table

It was found that locations with limestone and fissured and fractured and partly karst dolomite in high height and devoid of soil are associated with the chief karstic formations in western and northeastern parts of the area under study. Considering the ability of the mapping obtained in water resources management and determining the water balance components, including estimating penetration percentages (and consequently evaporation and groundwater storage) and identifying the location of the karstic resources compared to other existing methods much higher, the resulting base map, therefore, is used to locate and manage karstic water resources more economically and faster (Zagana et al. 2011).