Evaluation of Groundwater Vulnerability in the Upper Kelkit Valley (Northeastern Turkey) Using DRASTIC and AHP-DRASTICLu Models

: This study aimed to investigate groundwater vulnerability to pollution in the Upper Kelkit Valley (NE Turkey). For this purpose, vulnerability index maps were created using the generic DRAS-TIC and AHP-DRASTICLu models. The latter model was suggested by adding a parameter to the DRASTIC model and weighting its parameters with the analytical hierarchy process with the GIS technique. The results showed that areas with high and very high vulnerabilities are concentrated around the Kelkit Stream, which ﬂows from east to west in the central part of the study area. In contrast, areas with low and very low vulnerability classes are located in the northern and southern parts of the study area. To validate the model results, a physicochemical characterization of groundwater samples and their corresponding vulnerability index values were statistically compared using the Spearman correlation method. In addition, the single-parameter sensitivity method was applied to analyze the models’ sensitivities. Results revealed a stronger correlation between the vulnerability index values of the AHP-DRASTICLu model (compared to the DRASTIC model) in terms of sulfate (R 2 = 0.75) and chloride (R 2 = 0.76), while there was a slightly weaker correlation for the electrical conductivity (R 2 = 0.65) values of the groundwater samples. Sensitivity analysis indicated that the vadose zone, aquifer media, and land use are the most inﬂuential parameters responsible for the highest variation in the vulnerability index. Generally speaking, the results indicated that the AHP-DRASTICLu model performs better than the DRASTIC model for investigating groundwater vulnerability to pollution in the Upper Kelkit Valley.


Introduction
The ever-increasing population on a global scale has resulted in an increase in needs and therefore an increase in the resulting waste. This situation creates significant pressure on the environment, and this pressure is increasing day by day. In addition, climate change, in which the increasing population also has a large share, puts tremendous pressure on the environment [1]. This pressure mainly affects water resources [2]. Climate change projections reveal that the precipitation that Turkey receives in the future will decrease, leading to a decrease in water resources and, thus, the amount of usable water [3]. Model projections based on the IPCC-A2 pessimistic scenario show that there will be a 16% and 27% reduction in water potentials in Turkey by 2050 and 2075, respectively. This decrease in water potential due to the pressure of climate change [4] reveals the importance of using the existing water resources more effectively. In addition, keeping the quality of available water resources at a functional level is essential for this case.
Freshwater resources for human systems and ecosystems consist of groundwater and surface water resources. Groundwater supplies approximately half of the world's drinking water resources [5] and about 43% of all the water effectively used for irrigation [6]. Groundwater has a lower pollution potential compared to surface water, but polluted ISPRS Int. J. Geo-Inf. 2023, 12, x FOR PEER REVIEW 3 of 25 other large basins in Turkey, it is a region prone to water scarcity [3]. Regarding land use classes, more than half of the study area is of the farming activity class. The livelihoods of the majority of the local people depend on farming activity. Considering the population growth, an increase in the amount of agricultural activities, and therefore also in the chemical fertilizers used for product development, will be inevitable. Such activities on the land surface can disrupt essential natural ecosystem functions [46] and pose a threat to human health through a variety of exposure pathways [47]. In addition, groundwater is also used as drinking water in some parts of the study area. Therefore, studies on groundwater pollution in the Upper Kelkit Valley are a necessity for the proper management of groundwater and environmental planning. In addition, there are no studies on groundwater vulnerability in any terms in the region, and this study is a completely new attempt with the DRASTICLu model optimized with the combination of the AHP. Further, this new model, produced by adding land use criteria to the generic DRASTIC model and progressing the rates and weights of all criteria and subcriteria using the AHP, compared to previous models, will provide excellent convenience for groundwater vulnerability assessments in the areas similar to the UKV in terms of hydrogeoenvironmental conditions. Overall, this study aims to establish a groundwater vulnerability map of the Upper Kelkit Valley and provide proven bases for sustainable water management. In addition, modifying the generic DRASTIC model to create a more accurate groundwater vulnerability assessment model is also one of the aims of this study. Furthermore, this study aims to assess the relative importance of the parameters used in the models through sensitivity analysis for assessing aquifer vulnerability.

Description of the Study Area
The Upper Kelkit Valley (UKV) forms the upper parts of the Kelkit Valley, representing the upper parts of the Yeşilırmak Basin, one of Turkey's 25 largest river basins. The UKV extends between latitudes 40°03′45″-40°13′00″ N and longitudes 39°12′03″-39°42′13″ E and is located on the border of the Gümüşhane and Bayburt provinces (Figure 1). The topographic heights of the UKV, which covers an area of 445 km 2 , vary between 1357 and 2297 m a.m.s.l. (Figure 2a). In the study area with generally low land slope, the slopes vary between 0 and 18° in approximately 60% of its area ( Figure 2b, Table 1). The study area is located on the border of the Black Sea and the Eastern Anatolia Region and has transitional climatic conditions. Summers are cool; winters are quite cold and rainy. The meteorological measurements between 2003 and 2021 showed that the average annual precipitation in the Kelkit region was 300 mm and the average annual tem- The study area is located on the border of the Black Sea and the Eastern Anatolia Region and has transitional climatic conditions. Summers are cool; winters are quite cold and rainy. The meteorological measurements between 2003 and 2021 showed that the average annual precipitation in the Kelkit region was 300 mm and the average annual temperature was 6.5 • C [48]. The most common soil type in the study area is brown soil covering an area of 77% ( Figure 2c, Table 1). The rest consists of brown forest soil, bare rock, colluvium, noncalcic brown forest soil, alluvium, settlement, and floodplain ( Figure 2c). The UKV can be described as a rural area considering the settled population and land use classes. The UKV includes the Kelkit and Köse districts, which have populations of 44,068 and 7774, respectively [49]. The populations of these districts in 2007 were 41,664 (Kelkit) and 7270 (Köse), and their population growth rates are 12.64% (Kelkit) and 10.9% (Köse) [49]. This population increases even more during the summer seasons with the visits of relatives (i.e., grandparents, parents, siblings, children, and grandchildren) of citizens living in other cities or abroad. Arable land (35.55%) and agricultural lands (17.69%) are the land use classes that cover more than half of the UKV. Arable land, the most common land use type, is located around the Kelkit Stream, the main stream along the northeast-southwest direction of the UKV, in regions with an average slope of 4.  [50], and (d) geology-modified from [51].
Considering the temporal changes of the land use types in the study area [52][53][54][55], the 2000-2018 period witnessed a notable decrease in open spaces (29.96% and 18.61%) along with an increase in pastures (3.06% and 9.54%) and arable land (31.64% and 35.55%) (Figure 3, Table 2). The study area's residents mostly live from agriculture and animal husbandry activities. The increase in the agricultural and pasture areas in the study area  [50], and (d) geology-modified from [51].  The lithological basement unit in the study area is the Late-Carboniferous-aged Köse Granite [56,57], which is abundantly cracked and partially weathered [58] and crops out around the Köse settlement area in the north of the UKV. The Şenköy Formation spreading in different regions of the north and south of the study area [58], consisting of Liassic volcano sediments [59][60][61], overlies this unit unconformably. The Şenköy Formation is conformably overlain by the Berdiga Formation [60][61][62][63][64], which is Jurassic-Lower Cretaceous-aged shallow platform carbonates. The Berdiga Formation mainly consists of limestone, dolomite, and dolomitic limestone layers with conglomerate, sandstone, and silt successions that cover large areas in the northern and southern parts of the study area ( Figure 2d, Table 1). The Early-Middle-Eocene-aged Kelkit Formation (Figure 2d, Table  1), which outcrops in an east-west direction in the central parts of the study area, unconformably overlies the Şenköy and Berdiga Formations [51,65]. The formation mainly consists of volcanic and volcano-sedimentary rocks with abundantly fractured structures [58]. Alluviums (Figure 2d), which contain fragments of surrounding rocks from clay to coarse pebbles and generally outcrop along the Kelkit Stream, form the quaternary units in the study area [58].
According to the aquifer classification, considering the lithologies of the geologic formations in the UKV, the common aquifer types are the unconsolidated aquifer representing the alluvial material around Kelkit Stream, the semi-consolidated aquifer formed by the units of the Kelkit Formation surrounding this alluvial unit, and the karstic aquifer representing the Berdiga Formation scattered around the semi-consolidated aquifer unit in the northern and southern parts. In the study area, groundwater wells are mainly concentrated on unconsolidated and semi-consolidated aquifers around the Kelkit Stream. Groundwater abstracted from these aquifers and Kelkit Stream water is the primary water Considering the temporal changes of the land use types in the study area [52][53][54][55], the 2000-2018 period witnessed a notable decrease in open spaces (29.96% and 18.61%) along with an increase in pastures (3.06% and 9.54%) and arable land (31.64% and 35.55%) ( Figure 3, Table 2). The study area's residents mostly live from agriculture and animal husbandry activities. The increase in the agricultural and pasture areas in the study area through time can be explained by the increasing population growth rate and the fact that agriculture and animal husbandry are the main sources of income in the region . Irrigation water needed for the execution of farming activities is provided from the Kelkit Stream and its branches passing through the middle of the study area and from underground water wells. In addition, a large portion of the residents in the study area use the groundwater resources as drinking water. The lithological basement unit in the study area is the Late-Carboniferous-aged Köse Granite [56,57], which is abundantly cracked and partially weathered [58] and crops out around the Köse settlement area in the north of the UKV. TheŞenköy Formation spreading in different regions of the north and south of the study area [58], consisting of Liassic volcano sediments [59][60][61], overlies this unit unconformably. TheŞenköy Formation is conformably overlain by the Berdiga Formation [60][61][62][63][64], which is Jurassic-Lower Cretaceousaged shallow platform carbonates. The Berdiga Formation mainly consists of limestone, dolomite, and dolomitic limestone layers with conglomerate, sandstone, and silt successions that cover large areas in the northern and southern parts of the study area ( Figure 2d, Table 1). The Early-Middle-Eocene-aged Kelkit Formation ( Figure 2d, Table 1), which outcrops in an east-west direction in the central parts of the study area, unconformably overlies theŞenköy and Berdiga Formations [51,65]. The formation mainly consists of volcanic and volcano-sedimentary rocks with abundantly fractured structures [58]. Alluviums (Figure 2d), which contain fragments of surrounding rocks from clay to coarse pebbles and generally outcrop along the Kelkit Stream, form the quaternary units in the study area [58].
According to the aquifer classification, considering the lithologies of the geologic formations in the UKV, the common aquifer types are the unconsolidated aquifer representing the alluvial material around Kelkit Stream, the semi-consolidated aquifer formed by the units of the Kelkit Formation surrounding this alluvial unit, and the karstic aquifer representing the Berdiga Formation scattered around the semi-consolidated aquifer unit in the northern and southern parts. In the study area, groundwater wells are mainly concentrated on unconsolidated and semi-consolidated aquifers around the Kelkit Stream. Groundwater abstracted from these aquifers and Kelkit Stream water is the primary water source for agricultural and domestic use in the region. The Kelkit Stream, which flows through the study area in the northeast-southwest direction with a perennial annual average discharge rate of 0.403 m 3 /s [66], has an approximately 36 km main channel and its tributaries. Near the stream zone, the depth of groundwater is usually less than 1 m, and it rises to 99 m in areas with high altitudes (Figure 2a). Groundwater recharge to the aquifer is chiefly due to precipitation, subsurface inflows from the recharge area and Kelkit Stream channels, and agricultural irrigation.

Vulnerability Assessment Using DRASTIC and AHP-DRASTICLu Models
One of the most common methods used to estimate the vulnerability of aquifers to pollution is the US Environmental Protection Agency's DRASTIC [23] model. Although this model, which is a grid-cell-based overlay and index model, was not GIS-integrated at first, it was later modified by different researchers to various forms with the integration of GIS to further improve its pollution prediction capabilities in different physical environments [31][32][33]67]. In this study, both the DRASTIC model and the AHP-DRASTICLu model, which is a modified version of the DRASTIC model with the analytic hierarchy process, one of the multicriteria decisionmaking methods, were applied to the UKV to determine the groundwater vulnerability to pollution.
The DRASTIC method has two primary stages. First, the hydrogeological and geomorphological mappable parameters affecting groundwater pollution are classified. Then, using the determined coefficients, the parameters are weighted, and the vulnerability index is calculated [23].
Each letter of the word DRASTIC represents a hydrogeologically and morphologically mappable parameter [23]. These parameters are listed as follows: D: Depth to water table, R: Net recharge, A: Aquifer media, S: Soil media, T: Topography, I: Impact of vadose zone, C: Hydraulic conductivity. In the method, each parameter is rated (r), and weighted (w), and vulnerability maps are created on the GIS platform using the formula shown in Equation (1) determined for DRASTIC vulnerability index (DV i ) calculation [23].
In the rating (r) stage, each parameter of the DRASTIC model is considered a factor affecting pollution, and values ranging from 1 to 10 are assigned to these factors according to their relative pollution potential (Tables 3 and 4). "1" represents low potential, and "10" represents high potential. In the weighting (w) stage, weight values ranging from "1" (lowest importance) to "5" (highest importance) are assigned to each parameter that makes up the DRASTIC model, depending on their effects in transporting land-based pollutants to groundwater (Tables 3 and 4). Table 3. The pairwise comparison matrix and normalized weight (w i ) of the AHP-DRASTICLu model criteria and weight (w i ) of DRASTIC model criteria.  Table 4. The pairwise comparison matrix and normalized weight (w i ), standardized ratings (r i ) and total weight (r i × w i ) of the subcriteria of AHP-DRASTICLu and DRASTIC models.  Table 4. Cont.  In the AHP-DRASTICLu model, the land use classes parameter, an environmental factor that impacts groundwater pollution, was added to Aller's generic DRASTIC model. The rate and weight values of the main criteria (parameter) and their subcriteria were assigned with the analytical hierarchy process (AHP) method, which is one of the multicriteria decision-making methods. An AHP-DRASTICLu vulnerability index (AHP-DV i ) map was created by applying Equation (2) on the GIS platform to the rate (r) and weight (w) values calculated by the AHP for each main parameter of the AHP-DRASTICLu model. The methodology used in generic DRASTIC and AHP-DRASTICLu index maps is depicted in Figure 4.

Weight Assignment and Normalisation of DRASTICLu Criteria Using AHP
The analytical hierarchy process method [68], which is based on pairwise comparisons in line with expert opinions, provides great convenience in making complex decisions. The first stage of the AHP is determining the necessary criteria in line with the main target and forming pairwise comparison matrices (Equation (3)) of the selected criteria. Thus, the complex decision-making process between the criteria is reduced to a single level, and the relative importance values of the criteria are obtained. The degree of importance of the criteria against each other is assigned using values ranging from 1 to 9 (Table 5).
The most important subcriteria and their values are written in bold.
In the AHP-DRASTICLu model, the land use classes parameter, an environmental factor that impacts groundwater pollution, was added to Aller's generic DRASTIC model. The rate and weight values of the main criteria (parameter) and their subcriteria were assigned with the analytical hierarchy process (AHP) method, which is one of the multicriteria decision-making methods. An AHP-DRASTICLu vulnerability index (AHP-DVi) map was created by applying Equation (2)

Weight Assignment and Normalisation of DRASTICLu Criteria Using AHP
The analytical hierarchy process method [68], which is based on pairwise comparisons in line with expert opinions, provides great convenience in making complex decisions. The first stage of the AHP is determining the necessary criteria in line with the main target and forming pairwise comparison matrices (Equation (3)) of the selected criteria. Thus, the complex decision-making process between the criteria is reduced to a single  Table 5. Saaty's 1-9 scale of relative importance.

Scale
Judgment Explanation 1 Equally Two criteria contribute equally to the goal 3 Slightly Criterion 1 is slightly more important than criterion 2 5 Strongly Criterion 1 is strongly important compared to criterion 2 7 Very Strongly Criterion 1 is very strongly important compared to criterion 2 9 Extremely Criterion 1 is extremely important compared to criterion 2 2, 4, 6, 8 Intermediate Intermediate values between two adjacent numbers In Equation (3), a n shows the nth indicator unit, and a nn is the judgment matrix element. The second stage of the AHP methodology is to define the normalized weights using the geometric mean of the criteria (Equation (4)).
In Equation (4), W is the eigenvector and G m is the geometric mean of the ith row of the judgment.
The third and final stage of the AHP calculations is to test the consistency of the normalizations of the criteria. The consistency ratio (CR) (Equation (5)) is used in the consistency test. A CR less than 0.10 indicates that the normalized weight values of the criteria are consistent after pairwise comparisons, while a CR greater than 0.10 suggests that paired comparisons should be redesigned.
In Equation (5), CR is the consistency ratio, CI is the consistency index calculated using Equation (6), and RI is the random consistency index (see Table 6). In Equation (6), λ max is the maximum eigenvalue of the judgment matrix, which is calculated using Equation (7).
The main criteria (8 criteria) of the AHP-DRASTICLu model and pairwise comparison matrices of the subcriteria of these criteria, standardized rating (r), and normalized weight (w) values are given in Tables 3 and 4. As a result of the consistency tests of the pairwise comparisons of each criterion, it was determined that the CR values were lower than 0.10, and the values of the consistency test calculations (λ max , CI, and CR) are given in Table 7. Table 7. The number of criteria and their subcriteria (n), maximum eigenvalue of the judgment matrix (λ max ), consistency index (CI), random consistency index (RI), and consistency ratio (CR) for the criteria DRASTICLu model.

Data Preparation and Integration into a GIS Database
In this study, the basic datasets required to create the eight criteria layers used to determine the aquifer vulnerability to pollution have been obtained from different sources at various scales and resolutions (see Table 8). Each data layer was georeferenced within a GIS environment (ArcGIS 10.3.1) using the UTM projection system (Zone 37N) and WGS84 datum, and all operations (storage, integration, conversion, manipulation, analysis, and visualization) were performed using related plugins (such as 3D Analyst, Spatial Analyst, and Geostatistical Analyst) available in the software. The basic criteria layers in the GIS (D, R, A, S, T, I, C, Lu) used in this study consist of a grid format consisting of cells with a resolution of 25 × 25 m. The depth to water table refers to the distance from the surface to the groundwater table. The distance values measured in 40 different groundwater wells (Figure 1) scattered in the study area were used as basic data to create the depth to water table (D) layer in the UKV. Vector-based data transferred into the ArcGIS environment in point data format were interpolated using the kriging method, and the study area's depth to water table layer was created. Afterward, this grid layer was reclassified and divided into seven classes, and by assigning rating values (Table 4) to each class, "D" layers were created for the generic DRASTIC and AHP-DRASTICLu models (Figure 5a).  (Figure 1) scattered in the study area were used as basic data to create the depth to water table (D) layer in the UKV. Vector-based data transferred into the ArcGIS environment in point data format were interpolated using the kriging method, and the study area's depth to water table layer was created. Afterward, this grid layer was reclassified and divided into seven classes, and by assigning rating values (Table 4) to each class, "D" layers were created for the generic DRASTIC and AHP-DRASTICLu models (Figure 5a). The net recharge represents the annual total amount of water that infiltrates the vadose region and reaches the groundwater table [23]. The net recharge, one of the most The net recharge represents the annual total amount of water that infiltrates the vadose region and reaches the groundwater table [23]. The net recharge, one of the most important factors in transporting pollutants to groundwater, can be calculated by many different methods. In this study, the net recharge (R) layer of the study area was generated using the technique of Piscopo [70] (Equation (8)).
In Equation (8), RV is the potential recharge value, RA is rainfall amount, T is the topography (slope percentage), and SP is the soil permeability.
The precipitation amount layer of the study area was created by transferring the annual total average of the precipitation measured between 2003 and 2021 at Kelkit and Köse meteorological stations to the GIS environment. Depending on the precipitation amount of 300 mm/year throughout the study area, the RA value was assigned to the relevant criterion layer as "1" according to Table 9. The topography layer has been derived from the digital elevation model of the study area using the 3D analysis plugin in the ArcGIS environment. The digital elevation model of the study area was created by the "topo to raster conversion" plugin, which is available in ArcGIS, using vector data produced by digitizing the isohips of maps H42d, H43c, I42b, and I43a published by the Turkish Ministry of National Defense General Command of Maps. Afterward, this topography layer was divided into four different classes, as shown in Table 9. The T layer required for RV calculation was created by assigning the relevant rating value to the classes in the relevant range.
The soil permeability layers of the study area were obtained by digitizing the 1/25,000 scaled soil characteristics map of Turkey [50] in the ArcGIS environment. Depending on the texture characteristics of the soil types in the study area [71][72][73], they were divided into five classes: high, moderately high, moderate, low, and very low in terms of permeability. Then, the SP layer of the study area was created by assigning the rating values of each class (Table 9) to the relevant cells.
Finally, these RA, T, and SP layers were collected by the raster calculator plugin in ArcGIS, as shown in Equation (8). The last layer obtained was reclassified according to the RV ranges shown in Table 9. Then, "R" layers were created for the generic DRASTIC and AHP-DRASTICLu models by assigning the relevant rating values (Table 4) to the class ranges (Figure 5b).
The aquifer media represent the consolidated or unconsolidated material that forms the aquifer in the saturated zone [23]. The aquifer media (A) layer of the study area has been derived from the lithological units represented by the formations indicated in the 1/25,000-scaled geological map of the region [51]. The layer created in vector format with digitization in the ArcGIS environment was converted into the grid format with the help of the "feature to raster" plugin, and the "A" layers were created for the generic DRASTIC and AHP-DRASTICLu models by assigning the rating values (Table 4) to the relevant classes (Figure 5c).
The soil media (S) represent the uppermost layer of the vadose zone and are one of the most important control mechanisms in water infiltration into the ground and groundwater recharge [23,74]. The UKV's soil media (S) layer was derived from the 1/25,000-scaled soil characteristics map of Turkey [50] in vector format. Soil classes in the study area were classified according to the texture properties (clay loam for noncalcic brown forest and brown forest soils; loam for brown soil; sand for floodplain; sandy loam for alluvium and colluvium; thin or absent for settlement and bare rock) specified in some studies in the literature [71][72][73]. Afterward, this layer was converted into the grid format in the ArcGIS environment, and the soil media (S) layers of the study area were created for both models by assigning the rating values shown in Table 4 to the relevant classes (Figure 5d,e).
Topography refers to the slope of the land surface. The digital elevation model of the study area was used to create the topography (T) layers for the models of the UKV. Using the digital elevation model, the slope map was created with the help of the 3D Analysis tool in the ArcGIS software. Then, the slope map was divided into the classes indicated in Table 4, and the study area's topography (T) layers were created by assigning the rating values to the relevant classes (Figure 5f).
The vadose zone refers to the unsaturated or discontinuously saturated region between the groundwater table and the soil horizon [27]. The vadose zone material is one of the most significant mechanisms controlling transit time and the amount of pollutants and water leaking into the groundwater from the land surface [23]. The impact of the vadose zone (I) layers in the UKV has been derived using the vadose zone materials specified in the geological map of the study area [51]. The layer in vector format was converted into the grid format, then the "I" layers of the study area were created by assigning the rating values shown in Table 4 to the relevant classes (Figure 5g,h).
The hydraulic conductivity (C) is the capacity of the aquifer to transmit water [75]. The hydraulic conductivity layer of the UKV has been derived from lithological units [51,63] and well logs obtained from geotechnical studies [69] in the study area. The "C" layer of the study area was created by converting the layer in vector format into the grid format and assigning the rating values shown in Table 4 to the relevant classes (Figure 5i).
CORINE [55] land use land classes (LULCs) data in grid format were used as the basic data to create the land use classes (Lu) layer of the study area. The "Lu" layer of the study area was formed by assigning the rating values shown in Table 4 to the relevant classes ( Figure 5j).

Comparison and Validation of Generic DRASTIC and AHP-DRASTICLu Models
The similarity analysis of the generic DRASTIC and AHP-DRASTICLu model maps of the study area has been conducted by pairwise statistical comparisons under the Spearman correlation method [76]. In addition, the groundwater vulnerability indexes of both models have been validated with the Spearman correlation method using the sulfate, chloride, nitrate, and electrical conductivity of the groundwater samples collected from 20 wells ( Figure 6) in the study area in November 2022 and the vulnerability index values belonging to the well locations. scaled soil characteristics map of Turkey [50] in vector format. Soil classes in the study area were classified according to the texture properties (clay loam for noncalcic brown forest and brown forest soils; loam for brown soil; sand for floodplain; sandy loam for alluvium and colluvium; thin or absent for settlement and bare rock) specified in some studies in the literature [71][72][73]. Afterward, this layer was converted into the grid format in the ArcGIS environment, and the soil media (S) layers of the study area were created for both models by assigning the rating values shown in Table 4 to the relevant classes (Figure 5d,e). Topography refers to the slope of the land surface. The digital elevation model of the study area was used to create the topography (T) layers for the models of the UKV. Using the digital elevation model, the slope map was created with the help of the 3D Analysis tool in the ArcGIS software. Then, the slope map was divided into the classes indicated in Table 4, and the study area's topography (T) layers were created by assigning the rating values to the relevant classes (Figure 5f).
The vadose zone refers to the unsaturated or discontinuously saturated region between the groundwater table and the soil horizon [27]. The vadose zone material is one of the most significant mechanisms controlling transit time and the amount of pollutants and water leaking into the groundwater from the land surface [23]. The impact of the vadose zone (I) layers in the UKV has been derived using the vadose zone materials specified in the geological map of the study area [51]. The layer in vector format was converted into the grid format, then the "I" layers of the study area were created by assigning the rating values shown in Table 4 to the relevant classes (Figure 5g,h).
The hydraulic conductivity (C) is the capacity of the aquifer to transmit water [75]. The hydraulic conductivity layer of the UKV has been derived from lithological units [51,63] and well logs obtained from geotechnical studies [69] in the study area. The "C" layer of the study area was created by converting the layer in vector format into the grid format and assigning the rating values shown in Table 4 to the relevant classes (Figure 5i).
CORINE [55] land use land classes (LULCs) data in grid format were used as the basic data to create the land use classes (Lu) layer of the study area. The "Lu" layer of the study area was formed by assigning the rating values shown in Table 4 to the relevant classes (Figure 5j).

Comparison and Validation of Generic DRASTIC and AHP-DRASTICLu Models
The similarity analysis of the generic DRASTIC and AHP-DRASTICLu model maps of the study area has been conducted by pairwise statistical comparisons under the Spearman correlation method [76]. In addition, the groundwater vulnerability indexes of both models have been validated with the Spearman correlation method using the sulfate, chloride, nitrate, and electrical conductivity of the groundwater samples collected from 20 wells (Figure 6) in the study area in November 2022 and the vulnerability index values belonging to the well locations.

Sensitivity Analysis
An effective assessment tool is sensitivity analysis, which examines the contribution of individual variables and input parameters to the resultant output of an analytical model [77]. The sensitivity analysis validates and evaluates the uniformity of the analytical results of vulnerability maps [77]. It provides an efficient foresight to identify the most effective indicators of vulnerability and then explore how to deal with pollution and aquifer crisis management. In this study, the single-parameter sensitivity analysis, first introduced by Napolitane and Fabbri [77], was applied to check the accuracy of the vulnerability maps created by the DRASTIC and AHP-DRASTICLu models.
Single-parameter sensitivity analysis aims to compare the effective weights of different parameters with the theoretical weights assigned by the corresponding model [27]. The effective weight of each parameter was calculated by Equation (9) [77]: In Equation (9), W is the effective weighting of each parameter, P r and P w are the rating and weighting values of each parameter, respectively, and V is the overall vulnerability index.

Weighted Thematic Maps of the Criteria
The distances from the surface to the water table in the UKV vary between 0.5 m and 99.9 m. According to the generic DRASTIC method [23], the "Depth to water table" classes and covered areas in the UKV are 0.07% for 0-1.52 m, 0.42% for 1.52-4.57 m, 1.09% for 4.57-9.14 m, 9.14-4.  (Figure 5a). The depth to groundwater table indicates the distance at which the surface-borne pollutants will reach the groundwater, and groundwater is more vulnerable to pollution in areas with low distances [78]. Water and pollutants travel a longer pathway to reach the groundwater table in deeper groundwater table environments. This case allows more time for the pollutant to detoxify as it flows through the soil layers, indicating groundwater is less vulnerable to pollution [79]. In generic DRASTIC and AHP-DRASTICLu models, classes with low distances from the surface to groundwater have higher weight values, while classes with greater distance have lower weight values (Table 4).
Medium and high net recharge (R) classes are observed in the regions where the alluvial units formed by the Kelkit Stream accumulating the fragments that tore from the surrounding rocks around its bed, whereas the other parts generally have a very low net recharge, class "R" values ( Figure 5b). In both groundwater vulnerability models established for the UKV, low net recharge (R) classes have low weight values, whereas high net recharge classes have high weight values (Table 4). While 66.41% of the study area is represented by the very low net recharge (R) class, 29.24% is in the low class, 3.78% is in the medium class, and 0.57% is in the high class. The amount of recharge water leaking to the aquifer affects the amount of pollutants transported to the aquifer [80]. The higher the recharge water amount, the higher the vulnerability of groundwater to pollution and vice versa [81].
The geological formations in the study area were divided into "Alluvium", "Massive limestone", "Badded sandstone, limestone and shale", and "Igneous and Weathered metamorphic/igneous" classes according to the aquifer media (A) classification in the generic DRASTIC method in terms of their lithological properties. These classes cover 11.18%, 28.36%, 44.01%, 1.16%, and 15.30% of the study area, respectively.
The aquifer media classes in the UKV are as follows: the alluvial aquifer, by the alluvium units located around the Kelkit Stream; the "Badded sandstone, limestone, and shale" aquifer, by the conglomerate, sandstone, sandy limestone, shale, and tuff, which surround the alluvial aquifer and are referred to as the Kelkit Formation in the region; the "Massive limestone" aquifer, by the limestones of the Berdiga Formation, which is a karstic aquifer cropping out in topographically higher areas in the study area; the "Igneous" aquifer, by the granites covering a relatively small area in the northeast of the study area; and the "Weathered metamorphic/igneous" aquifer, by the "Şenköy Formation" units, which cover areas in the northern and southern parts of the study area, which generally consist of basaltic and pyroclastic rocks with abundant cracks and fractures (Figures 2d and 5c). While the alluvium aquifer class has the highest weight value in the aquifer media criterion for both models, the lowest value belongs to the "igneous" aquifer class ( Table 4). The flow path of water below the surface is controlled by factors such as the hydraulic conductivity of the aquifer and gradient [82]. Since the fine-grained and low-permeability aquifer material is more resistant to the infiltration of pollution into the aquifer, the pollution vulnerability is lower in these systems.
The most common soil class in the UKV is brown soil with a loam texture, which covers approximately 3/4 of the study area. This class, which has a rating of "5" according to the generic DRASTIC method, has a rating of "0.047" according to the AHP-DRASTICLu model ( Table 4). The rest of the soil classes with areas they cover in the UKV are noncalcic brown forest and brown forest soil with a clay loam texture of 6%, floodplain with a sand texture of 0.94%, alluvium and colluvium with a sandy loam texture of 7.72%, and bare rock and settlement represented by thin or absent classes of 6.21%. In soils with more coarsegrained materials, such as sand and gravel, the water infiltration capacity will be higher, and the vulnerability of groundwater to pollution is higher in these systems [83]. However, the presence of organic materials, such as clay and organic matter, creates resistance against water infiltration to the ground, reducing pollution vulnerability in groundwater [33,84]. Among the soil classes in the UKV, the highest vulnerability rating value belongs to the bare rock and settlement classes, represented by "thin or absent" for both models. In contrast, the lowest rating values are for noncalcic brown forest and brown forest soils depicting the clay loam texture spread in different regions of the UKV (Table 4 and Figure 5d,e).
The topography layer of the UKV was classified into five groups according to land slope (%) ranges: 0-2%, 2-6%, 6-12%, 12-18%, and >18%. These classes occupy the areas of 6.28%, 18.61%, 21.08%, 15.02%, and 39.01% in the UKV, respectively (Table 1). The topography criterion is one of the most critical factors controlling the water residence time on the land surface [85]. In the land with a low slope, the water stays on the surface longer, so the amount of surface water infiltrating underground will be greater. Thus, the amount of pollutants that can leak into the groundwater also increases [81]. In the northern and southern parts of the study area, the slope of the land is high in topographically high regions, and the rating values of the topography classes are low in both generic DRASTIC and AHP-DRASTICLu models in these areas (Figure 5f). In the middle parts of the study area, the land slopes along the east-west extension belt are lower, and the topography class rating values are higher for both models in these areas (Figure 5f).
The vadose zone materials in the UKV were grouped into five different classes, namely "Alluvium", "Massive limestone", "Badded limestone, sandstone and shale", "Igneous", and "Weathered Metamorphic/igneous", as specified in the DRASTIC method. These classes cover 11.18%, 28.36%, 44.01%, 1.15%, and 15.30% of the UKV, respectively. The fact that the vadose zone consists of units with a high permeability, such as sand and gravel, allows pollutants to reach the groundwater in larger quantities and in a shorter time. This situation indicates that the vulnerability of groundwater to pollution in these regions is higher [23]. Among the vadose zone classes in the UKV, the most vulnerable class for groundwater pollution is "Alluvium", and the rating values are "8" for the generic DRASTIC model and "0.438" for the AHP-DRASTICLu model. In contrast, the class with the lowest vulnerability is "Igneous", and the rating value of this class is "4" for the generic DRASTIC model and "0.047" for the AHP-DRASTICLu model (Figure 5g,h).
The movement of water and pollutants through the aquifer is greatly affected by hydraulic conductivity. The higher the hydraulic conductivity of an aquifer, the more pollutants could be transmitted to groundwater and, therefore, the higher the groundwater vulnerability to pollution [10,33,85]. The hydraulic conductivity values of the UKV were classified into three different classes: 12.225-28.525 m/d, 28.525-40.75 m/d, and 40.75-81.5 m/d. While the highest hydraulic conductivity class has the highest rating values of vulnerability for both the generic DRASTIC (8) and AHP-DRASTICLu (0.243) models, the class with the lowest hydraulic conductivity values has the lowest ratings for both models (4 for the DRASTIC and 0.083 for the AHP-DRASTICLu) ( Table 4). The areas with the highest hydraulic conductivity rating values and, therefore, the areas most vulnerable to pollution in terms of hydraulic conductivity, are located around the Kelkit Stream, which is the main stream of the study area (Figure 5i). This class covers 11.17% of the study area. However, the lowest vulnerable areas in terms of hydraulic conductivity cover approximately 3/4 of the UKV.
Land use classes can be a source of pollutants and can control the way water [86] and pollutants reach the groundwater [10,87,88]. Anthropogenic effects, such as domestic and industrial wastes, improper design of septic tanks, and agricultural activities, especially pesticides, may affect the groundwater quality and cause it to become polluted [38]. A large part of the UKV consists of farming activity classes (35.5% of UKV is arable land, and 17.7% of UKV is agricultural area), and these classes are located at the central part of the study area along an east-west extending belt with a low land slope (Figure 3d). These classes are highly vulnerable to pollution, and the rating values in the AHP-DRASTICLu model are "0.192" for "agricultural area" and "0.135" for "arable land" ( Table 4). The highest rating value of the land use classes criteria belongs to the settlement with "0.257". The rating values for the other land use classes, "forestry", "mine, dump site", "natural grassland", "open field", "pasture", "shrub", and "water body", are "0.064", "0.126", "0.046", "0.053", "0.055", "0.054", and "0.018", respectively ( Table 4). The spatial distributions of rated land use classes are shown in Figure 5j.

The Groundwater Vulnerability Maps of DRASTIC and AHP-DRASTICLu Models
"Depth to water table", "Net recharge", and "Impact of vadose zone" criteria, compared to the other criteria used in the DRASTIC and AHP-DRASTICLu models, in addition to the "Land use" criterion with its impact in AHP-DRASTICLu model, are more decisive factors on groundwater vulnerability to pollution according to their theoretical weights. Many previous studies have shown that these criteria are dominant factors in groundwater vulnerability to pollution [10,31,33,88]. The generic DRASTIC index (DV i ) values of the study area obtained as a result of the analysis range between "64" and "193", and the AHP-DRASTICLu index (AHP-DV i ) values vary between "0.042" and "0.302" (Figure 6). These vulnerability index maps of the UKV were reclassified into five categories (very high vulnerability, high vulnerability, moderate vulnerability, low vulnerability, and very low vulnerability) to help better understand the spatial distribution patterns of the different vulnerability categories. The areal distributions (% area of the UKV) of these classes were 1.42%, 5.48%, 7.80%, 34.31%, and 50.98%, respectively, according to the generic DRASTIC vulnerability index classification, and they were 3.90%, 8.02%, 31.78%, 29.86%, and 26.44%, respectively, according to the AHP-DRASTICLu vulnerability index classification.
As a result of both models, classes with high to very high vulnerabilities were observed along the east-west belt formed along the Kelkit Stream in the study area ( Figure 6). The areas represented by these classes are generally land use classes related to farming activity with a very low topographic slope (with averages of 2 • for the DRASTIC model and 4 • for AHP-DRASTICLu model), dominated by alluvium and floodplain soil classes and areas with relatively shallow water table depths where residential areas are the dominant land use classes. Similar results were reported in the previous studies on groundwater vulnerability in northern Turkey [85] and in the arid/semi-arid regions with climatic conditions similar to the UKV [33,67,89].
The moderate vulnerability class of the AHP-DRASTICLu model covers a larger area than the DRASTIC model. These classes mainly extend to the high vulnerability class borders with the east-west trending hills in the northern and southern parts of the study area, where the topographical slope begins to increase relatively. The average slope of the moderate vulnerability class of the DRASTIC model is 7 • , while it is 6 • in the AHP-DRASTICLu model. The dominant aquifer type in these areas consists of clastic rocks containing fine-grained materials. The predominant soil class is "brown soil", and the most common land use classes are "arable land" and "agricultural area" in the moderate vulnerability areas of the UKV.
In the northern and southern parts of the study area, the vulnerability varies from low to very low on the east-west trending hill range, where the topographic heights and slopes are relatively high ( Figure 6). In areas representing this class, the distances to the groundwater table are relatively greater (an average of 34 m for the DRASTIC model, and an average of 38 m for the AHP-DRASTICLu model).

Models' Comparison and Validation
In order to test the accuracy of the vulnerability index maps obtained, 20 different groundwater samples scattered over the area in November 2022 were gathered (Figure 6), and their physicochemical contents were analyzed. The analysis results and vulnerability index values were compared using the Spearman correlation method.
The sulfate concentrations of the groundwater samples collected from the study area varied from 15.12 mg/L to 115.14 mg/L, the chloride concentrations of the samples ranged between 4.89 mg/L and 58.63 mg/L, the nitrate concentrations of the samples ranged between 0.51 mg/L and 200.45 mg/L, and the electrical conductivity values ranged between 499 µS/cm and 1450 µS/cm. While calcium sulfates, such as gypsum and anhydrite, are the most common sources of sulfates in natural waters, barium sulfates, strontium sulfates, the oxidation of pyrite minerals, volcanoes, geothermal fields, the burning of fossil fuels, and the oxidation of sulfates in soil or organic wastes are also sources of sulfate in natural water [90,91]. Apart from these factors, the sulfate concentration in surface water and groundwater may increase due to industrial and domestic wastes and sulfatecontaining fertilizers used in agriculture [92][93][94]. Chloride leaching from chemical fertilizers in agricultural soils or wastewater discharged to the land surface may increase the chloride concentration in groundwater [95]. The most common nitrate sources in groundwater are inorganic and organic fertilizers used in agricultural activities, organic materials in the soil, industrial wastes, solid waste disposal areas, highways, and sewage waters [96][97][98]. The electrical conductivity of water represents the dissolved ion concentration in the water. The dissolved ion concentration in water with high electrical conductivity is high and vice versa [90,99].
As a result of the correlation analysis, it was observed that there was a strong correlation between the vulnerability index values of both models and the sulfate, chloride, and electrical conductivity contents of the groundwater samples, whereas the nitrate values and vulnerability index values were not correlated. In the UKV, the correlation coefficients between sulfate concentrations of groundwater samples and generic DRASTIC and AHP-DRASTICLu index values are "0.661" and "0.752", respectively; the correlation coefficients between chloride concentrations and generic DRASTIC and AHP-DRASTICLu index values are "0.684" and "0.758", respectively; the correlation coefficients between nitrate concentrations and generic DRASTIC and AHP-DRASTICLu index values are "−0.245" and "−0.151", respectively; and the correlation coefficients between the electrical conductivity and DRASTIC and AHP-DRASTICLu index values are calculated as "0.678", and "0.652", respectively (Table 10). The correlation coefficient between the DRASTICLu index values and the generic DRASTIC index values is "0.931". Considering the correlation analysis results, it was seen that the AHP-DRASTICLu model provides more accurate results than the DRASTIC model in terms of the sulfate and chloride contents of the groundwater except for the electrical conductivity in the UKV. In many studies conducted to determine the groundwater vulnerability to pollution, the nitrate [31,100,101], sulfate [102][103][104], and chloride [29,105] concentrations of groundwater were used to validate the model results. It was seen that the vulnerability index values were highly correlated with sulfate in some areas [106,107], chloride in some areas [35,108], and nitrate in some areas [36,109]. These differences suggest groundwater vulnerability models respond differently in different climatic and hydrogeoenvironmental areas.

Single-Parameter Sensitivity Analysis
Single-parameter sensitivity analysis compares "theoretical" and "effective" weights. Here, the effective weight is a function of the value of the single parameter compared to other parameters and the weight the model assigns [27]. The effective weights of the parameters of the DRASTIC and AHP-DRASTICLu models have deviated from theoretical weights (Table 11). In the DRASTIC and AHP-DRASTICLu models, the impact of the vadose zone parameter has a higher effective weight value (30.11% and 40.64%, respectively) than the other parameter, as well as a higher value than the relevant theoretical weights. It is seen that the impact of the vadose zone parameter is the most influential parameter for the sensitivity assessment. Further, the effective weights of aquifer media (10.83%) and hydraulic conductivity (15.34%) parameters in the DRASTIC model and the aquifer media parameter (13.74%) in the AHP-DRASTICLu model show higher values than their theoretical weights (Table 11). The effective weights of all other layers are lower than their theoretical weights. These results show that the impact of the vadose zone parameter has a very important role in both models, and it is important to obtain accurate and detailed information on this factor. In addition, the high effective weight value (14.55%) of the land use classes parameter in the AHP-DRASTICLu model has shown that including this parameter in the model is quite accurate.

Conclusions
In the present study, the groundwater vulnerability maps of the Upper Kelkit Valley (UKV) have been derived using the generic DRASTIC model and the AHP-DRASTICLu model, which is a modified model that adds the land use class criterion and the weighting of the criteria by the analytical hierarchy process method. According to the vulnerability index values obtained from both models, the study area was reclassified into five classes: very low, low, moderate, high, and very high vulnerability. Considering the spatial distributions of the vulnerability index, the areas covered by the very high, high, and moderate vulnerability classes in the AHP-DRASTICLu model index map in the UKV are larger than in the generic DRASTIC model index map. Both models show that the high to very high vulnerability classes are extended along a belt around the Kelkit Stream in the middle parts of the study area. In contrast, the lower vulnerability classes fall in the areas where the distance to the groundwater is relatively higher, along the east-west trending hills in the north and south of the study area. Particularly, the high and very high vulnerability classes are the areas around the Kelkit district, which is the largest settlement of the study area and where agricultural activities are intense . Considering these factors create a high pollution potential for groundwater, it is recommended that water and land use management planning should be more careful with an integrated management plan by all related persons, institutions, and groups, from water users to decision makers.
In order to validate the model results, the physicochemical contents of the water samples taken from different groundwater wells scattered in the study area and the vulnerability index values of each well location were statistically compared using the Spearman correlation method. The comparison results showed a strong correlation between the groundwater samples' electrical conductivity, sulfate, and chloride concentrations and the vulnerability indexes of both models. In addition, the correlation results showed that the AHP-DRASTICLu model is more accurate than the generic DRASTIC model for groundwater vulnerability assessment of the UKV. Furthermore, the sensitivity analysis revealed that the impacts of the vadose zone, aquifer media, and land use are the most influential parameters responsible for the highest variation in the vulnerability index. These parameters are significant factors as they have a principal role in determining accurate and detailed information for water and environmental planning.
The AHP method reduces the bias due to standard weights in the generic DRASTIC model when combined with the DRASTIC model, one of the most effective methods for vulnerability estimation in groundwater. In addition, adding the land use criterion to the model simulates better the pressures that may cause groundwater pollution. In this sense, with the AHP-DRASTICLU model used in this study, the groundwater vulnerability map for the study area was simulated to be closer to the actual pollution vulnerability.
The fact that the AHP-DRASTICLu model, which was established by modifying the generic DRASTIC model, provides more accurate results indicated that the AHP-DRASTICLu model could be used more efficiently instead of the generic DRASTIC model. However, the AHP-DRASTICLU model may give different results in different climatic and hydrogeoenvironmental areas. Therefore, AHP analyses should be rearranged according to regional conditions. Funding: This work was supported by the Research Fund of the Bayburt University, Project Number: 2022/69002-13.
Data Availability Statement: The hydrological and hydrogeologic input data used for the DRASTIC and AHP-DRASTICLu models presented in this study are available upon request from the author.