Groundwater Vulnerability Assessment—Case Study: Tirana–Ishmi Aquifer, Albania

: This paper discusses the groundwater vulnerability to pollution assessment for the Tirana–Ishmi alluvium aquifer, Albania. Economic activities, municipal wastewater discharged into rivers and groundwater overexploitation threaten to pollute the groundwater. Based on the aquifer characteristics and the available data, SINTACS was selected as the most realistic assessment model. The SINTACS parameters’ rates assigned to the aquifer’s characteristics (water table depth, infiltration, unsaturated zone, soil media, aquifer media, hydraulic conductivity, topography) were adapted to the local features, followed by GIS vulnerability mapping. Statistical analysis indicates that the unsaturated zone, hydraulic conductivity and aquifer media have the highest influence on groundwater vulnerability, whereas topography has the lowest influence. Validation through sensitivity analysis and nitrates content confirms the rational selection of the SINTACS model and the reliability of the study’s outputs. The most vulnerable areas to pollution are the recharge zones, followed by the highly urbanized Tirana City area, characterized by high levels of groundwater extraction rate and wastewater discharged into the rivers. The paper, being the first completed groundwater vulnerability assessment of the study area, could serve as a basis for a scientific–based groundwater management that should be considered in local territory planning.


Introduction
The present paper discusses the groundwater vulnerability (GWV) assessment of the Tirana-Ishmi aquifer, in Albania, and the related mapping using a Geographical Information System (GIS).The concept of aquifer vulnerability expresses the aquifer's susceptibility to be affected by contaminant substances from the land surface [1][2][3].This susceptibility is defined in the context of the standard hazard-pathway-receptor model [4,5].The vulnerability is a relative, non-measurable and dimensionless aquifer attribute indicating the aquifer's sections where contamination is most likely to occur [6,7].The GWV can be intrinsic or specific.The intrinsic vulnerability depends on the hydrogeological characteristics of an aquifer, while the specific vulnerability is based on the risk of the contamination generated by anthropogenic activities [7].Based on the concept of GWV and the hydrogeological factors influencing it, the first vulnerability maps were prepared at the end of the 1960s in countries such as the USA, Czechoslovakia, etc. [8], and in 1970 in France [9].Aiming at groundwater protection and its optimal usage, the hydrogeologists' efforts worldwide were focused on the detailed study of GWV and its influencing factors [10], as well as in the preparation of GWV maps, which should be taken into account by decision makers for the groundwater's protection from pollution and territory planning [11,12].
Based on the characteristics of the study area and the available data and information, the intrinsic vulnerability is the appropriate choice for assessing the GWV of the Tirana-Ishmi aquifer.Amongst numerous intrinsic methods [13], SINTACS [14][15][16] was selected as the most realistic one for the purposes of this study.The parameters taken into consideration by SINTACS are water table depth, infiltration, unsaturated zone, soil media, aquifer media, hydraulic conductivity and terrain slope.The assigned rates to these parameters were adapted to local features.
The GWV was mapped using a GIS [16,17].The used data were submitted to a statistical analysis [18], which indicates that the unsaturated zone, hydraulic conductivity and aquifer media have the highest influence on the GWV, whereas the topography has the lowest influence.Validation through sensitivity analysis (SA) [19] and nitrates content [20] confirmed the rational selection of the used GWV assessment method and the reliability of the study's outputs and outcomes.The vulnerability map, which is the main output of this study, indicates that the areas most vulnerable to pollution are the recharge zones, followed by the highly urbanized Tirana City area, characterized by high levels of both groundwater extraction rate and wastewater discharged into the rivers.
The study area is located in Middle Albania [21].It represents a flat alluvial valley bordered by an open flat terrain northwestward and by molassic hills from the other sides (Figure 1) [21,22].The aquifer is composed of Quaternary alluvial deposits that outcrop only in the rivers' beds in the eastern edge of the study area [21,23,24].The area is home to nearly 650,000 inhabitants, who use the groundwater for drinking, domestic, industrial and irrigation purposes.The local population in the southern part of the area has significantly grown in the last 30 years due to uncontrolled demographic movements.At the same time, the extracted groundwater amount in this area has substantially intensified.[25].Municipal wastewater is discharged into the middle and lower course of the local rivers [25,26].An exception is made for some villages in the central and northern part of the lowland, whose sewage is discharged into private septic holes or tanks.As of the time of preparation of this paper, there is no wastewater treatment plant (WWTP) in the study area [25], although the construction of a WWTP for Tirana City is already planned [27].The wastewater flowing into the rivers poses a threat to the population's health, due to the risk of groundwater pollution [28][29][30], which is associated with the hydrogeological and hydrological features of the study area, as well as with the location and characteristics of the groundwater pollution sources.
The GWV assessment in Albania has been subject to four scientific papers, three of which discuss the GWV in other regions of the country [31][32][33], while the fourth paper aims rather to show how to use GIS technology for groundwater protection from pollution [33].None of these papers use either a statistical analysis or a validation of the GWV assessment.The fourth paper [34], whose study area includes the Tirana-Ishmi aquifer, has limitations regarding the GWV assessment as it does not provide the hydrogeological settings, does not indicate the potential sources of groundwater pollution, does not present the map of each DRASTIC parameter, does not validate the used assessment, etc.Therefore, the present article can be considered as the first completed scientific paper attempting to integrally assess the vulnerability of the Tirana-Ishmi aquifer.
Based on the reliability of its outputs and expected outcomes, this article could serve as a reliable basis for preparing a scientific-based management plan of the aquifer and for inclusion of the GWV assessment in the development plans of the study area [16,17].The absence of a GWV discussion and elaboration in the recent "Management plan of the Ishmi River basin for the period 2024-2029" [25], approved by the Albanian government [35], indicates that the decision makers are not yet aware either of the importance of GWV in the river basins' management plans and other spatial development plans, or of a closer involvement of hydrogeology experts in the preparation of these plans.Based on the reliability of its outputs and expected outcomes, this article could serve as a reliable basis for preparing a scientific-based management plan of the aquifer and for inclusion of the GWV assessment in the development plans of the study area [16,17].The absence of a GWV discussion and elaboration in the recent "Management plan of the Ishmi River basin for the period 2024-2029" [25], approved by the Albanian government

Study Area 2.1. General Features
The Tirana-Ishmi aquifer lies in a flat valley oriented in the SE-NW direction (Figure 1).Its surface area is roughly 190 km 2 and the length is about 37 km, while the width varies roughly from 3 km to 10 km [36].The approximate altitude above the sea level (a.s.l) is 100 m in the south, 30 m in the center and 2.0 m in the north.The altitude decreases from east to west, by roughly 50 m a.s.l. in the south, 20 m a.s.l. in the center and 11-12 m a.s.l. in the north.The height of the hilly ranges surrounding the study area varies from 310 to 390 m a.s.l.Eastward of the hills that border the study area from the east lies the Kruje-Dajti carbonate mountain range, whose maximum height is 1613 m a.s.l.[22,36].The small rivers of Lana, Tirana, Terkuza, Zeza and Droja, all of which originate at the Kruje-Dajti mountain range, traverse the study area and form the Ishmi River, which is a typical lowlands river.There are no water springs in the study area.The closest springs are located in the Kruje-Dajti karst massif [21,23,24,37,38].Soft and wet winters and hot and dry summers characterize the climate in the study area [39].The annual average, the minimum absolute and maximum absolute temperatures are 15.0 • C, −12.9 • C and 41.5 • C, respectively.The annual average rainfall is roughly 1250 mm/year [39][40][41].
The main urban centers in the study area include Tirana City and the Kamza, Fushe Kruje and Mamurras towns (Figure 1).A considerable number of population live in agricultural areas.The population of Tirana and Kamza in the last 30 years has drastically grown due to uncontrolled population movements.The population in the southern half of the study area is about 600,000 inhabitants [25].The industrial facilities are mainly located in the southwest of the area.They comprise car reparation units, facilities for the assemblage of laundry and dishwasher detergents, facilities for the production of alcoholic and non-alcoholic beverages, wood, milk, cereals and meat processing plants, etc.The wastewater of these facilities is discharged into the Lana and Tirana Rivers [25].The sewage flows into the same sewers as rain and then drains into the nearby rivers by gravity, thanks to the favorable terrain slope [25].The sewerage flows into the middle and lower course of the closest rivers.Currently, the urban wastewater of the Tirana-Kamza area flows into the Lana and Tirana Rivers, that of Fushe Kruje flows into the Zeza River and that of Mamurras into the Droja River [25,42].The urban solid waste of Tirana City and its neighboring area is disposed of into a landfill and does not pollute either the surface waters or the groundwater.The use of fertilizers in agriculture over the last 30 years is almost absent.
The extracted groundwater amount from hydrogeological wells is about 970-1000 L/s for drinking purposes, 200-300 L/s for technological purposes and about 100-200 L/s for irrigation.The amount extracted from private wells is not known.The total amount varies from 1270 to 1500 L/s [25].The annual groundwater utilization rate in the southern part of the study area varies from 0.85 to 0.95 [43].The aquifer supplies drinking water to almost all the villages of the study area and a part of Tirana City.The wells in the territory of this city supply roughly 300-490 L/s [25].The remaining water amounts needed for Tirana are supplied by other sources located outside the study area [43].

Geological Settings
The study area is located in the Tirana-Ishmi syncline [23,24,44].The hills surrounding the study area from the east, west and south consist of Upper Miocene molasses composed of Serravalian (N 2  1 s), Tortonian (N 3 1 t) and Messinian (N 3 1 m) age deposits (Figure 2) [22] represented by intercalations of coarse-grain sandstone layers with siltstone and clay layers [22,45,46].The upper geological section of the area consists of prolluvial (Qp-h) and alluvial (Qh) Quaternary deposits.The thickness of the prolluvial deposits varies from 3.0 m to 8.0 m [47].They are mainly composed of gravel, sand and siltstone originating from the erosion of the hilly slopes and river and stream beds.They extend mainly to the east and south of the study area.In the western part of the area, they form a narrow strip along the foothills.[22] (Figure 2).The alluvial deposits (Qh), which cover the central and western parts of the study area, are represented by the terraces of the rivers that cross the area.The alluvial deposits are composed of gravels, sands and clays [21,45].
m to 8.0 m [47].They are mainly composed of gravel, sand and siltstone originating from the erosion of the hilly slopes and river and stream beds.They extend mainly to the east and south of the study area.In the western part of the area, they form a narrow strip along the foothills.[22] (Figure 2).The alluvial deposits (Qh), which cover the central and western parts of the study area, are represented by the terraces of the rivers that cross the area.The alluvial deposits are composed of gravels, sands and clays [21,45].The thickness of the alluvial deposits increases from south to north and from east to west.In the south to north direction, this thickness varies from nearly 5 m in Tirana City to 100-120 m in the northwestern edge of the study area (Figure 3).Meanwhile, in the east to west direction, in the Fushe Kruje-Mamurras sector, this thickness varies from 20 m in the east to 100-120 m in the west [47].The thickness of the alluvial deposits increases from south to north and from east to west.In the south to north direction, this thickness varies from nearly 5 m in Tirana City to 100-120 m in the northwestern edge of the study area (Figure 3).Meanwhile, in the east to west direction, in the Fushe Kruje-Mamurras sector, this thickness varies from 20 m in the east to 100-120 m in the west [47].

Hydrogeological Settings
Several scientific papers describe in detail the hydrogeological features of the aquifer under study [48][49][50][51][52][53][54][55][56][57][58], which consists of a single gravelly layer, with the exception of some sectors where the gravels form two or more layers separated by clayey lentils (Figure 3) [54,55].Based on the lithological data from 35 wells [50,54,57,58], we applied Groundwater Modeling System (GMS) software [59] to create a 3D model of the aquifer (Figure 4).The gravelly layers have good hydraulic connections between them.Their thickness increases from south to north.In the west of Tirana City, this thickness is about 5 m, while in the north of the aquifer, it reaches 60-70 m.The only exception consists of a limited area in the southwest of Tirana, where the thickness of the water-bearing gravels reaches 31 m.The gravelly aquifer is overlaid by a clayey layer, the thickness of which increases from 2-3 m in the south to 20 m to 30 m in the north of the study area [54].The gravelly deposits outcrop on the surface only in the recharge zones, which are located in the eastern part of the aquifer [49,54,56,58] (Figure 4).The aquifer is unconfined in the southern sector (Tirana City and its neighboring area) and confined in the central and northern sectors.In the northwest of the study area, the water wells are artesian.In the sectors where the aquifer is unconfined (or phreatic), it is recharged by the infiltration of atmospheric precipitation and the Lana, Tirana, Terkuza, Zeza and Droja Rivers.The possibility of the aquifer's recharging from the underlying mollassic deposits is not excluded.Although these deposits have a low water-bearing capacity, their water could drain through sandstone fractures due to high pressure [57].The main groundwater flow direction in the aquifer is southeast-northwest [21], with a hydraulic slope gradient of I = 0.007-0.009[57].
Information on the aquifer's hydraulic properties is based on the Albanian Geological Survey (AGS) studies and other scientific papers [50,57,58], which are all based on the 130 wells' pumping data collected by the AGS during the period 1964-2011 [57,58].The hydraulic conductivity in the recharge zones of the Tirana and Terkuza Rivers varies from 180 to 300 m/day [57,58].Similarly, the recharge zones of the Zeza and Droja Rivers have a high hydraulic conductivity, too [60].In the southwestern edge of Tirana City (Selita Vogel) and in the other main wellfields (Laknas-Berxulle, Rinas-Valias and Larushk-Fushe Kruje), located in the center of the study area, the hydraulic conductivity varies from 100 to 300 m/day.Around these wellfields, the hydraulic conductivity varies from 50 to 90 m/day, whereas in the other parts of the aquifer, it is lower than 50 m/day [57,58].The

Hydrogeological Settings
Several scientific papers describe in detail the hydrogeological features of the aquifer under study [48][49][50][51][52][53][54][55][56][57][58], which consists of a single gravelly layer, with the exception of some sectors where the gravels form two or more layers separated by clayey lentils (Figure 3) [54,55].Based on the lithological data from 35 wells [50,54,57,58], we applied Groundwater Modeling System (GMS) software [59] to create a 3D model of the aquifer (Figure 4).The gravelly layers have good hydraulic connections between them.Their thickness increases from south to north.In the west of Tirana City, this thickness is about 5 m, while in the north of the aquifer, it reaches 60-70 m.The only exception consists of a limited area in the southwest of Tirana, where the thickness of the water-bearing gravels reaches 31 m.The gravelly aquifer is overlaid by a clayey layer, the thickness of which increases from 2-3 m in the south to 20 m to 30 m in the north of the study area [54].The gravelly deposits outcrop on the surface only in the recharge zones, which are located in the eastern part of the aquifer [49,54,56,58] (Figure 4).The aquifer is unconfined in the southern sector (Tirana City and its neighboring area) and confined in the central and northern sectors.In the northwest of the study area, the water wells are artesian.In the sectors where the aquifer is unconfined (or phreatic), it is recharged by the infiltration of atmospheric precipitation and the Lana, Tirana, Terkuza, Zeza and Droja Rivers.The possibility of the aquifer's recharging from the underlying mollassic deposits is not excluded.Although these deposits have a low water-bearing capacity, their water could drain through sandstone fractures due to high pressure [57].The main groundwater flow direction in the aquifer is southeast-northwest [21], with a hydraulic slope gradient of I = 0.007-0.009[57].
Information on the aquifer's hydraulic properties is based on the Albanian Geological Survey (AGS) studies and other scientific papers [50,57,58], which are all based on the 130 wells' pumping data collected by the AGS during the period 1964-2011 [57,58].The hydraulic conductivity in the recharge zones of the Tirana and Terkuza Rivers varies from 180 to 300 m/day [57,58].Similarly, the recharge zones of the Zeza and Droja Rivers have a high hydraulic conductivity, too [60].In the southwestern edge of Tirana City (Selita Vogel) and in the other main wellfields (Laknas-Berxulle, Rinas-Valias and Larushk-Fushe Kruje), located in the center of the study area, the hydraulic conductivity varies from 100 to 300 m/day.Around these wellfields, the hydraulic conductivity varies from 50 to 90 m/day, whereas in the other parts of the aquifer, it is lower than 50 m/day [57,58].The average specific discharge within these wellfields varies from 10 to 15 L/s/m, while the maximum specific discharge in some wells varies from 20 to 30 L/s/m [58].
net recharge and the molassic deposits is roughly 100 L/s [63].The total inflow (710 L/s) is lower than the total outflow, which varies from 1270 to 1500 L/s [25].The increase in the groundwater outflow compared to the natural inflow occurs thanks to induced filtration due to the groundwater recharging directly from the rivers and the lowering of the dynamic groundwater level.
The groundwater of the Tirana-Ishmi aquifer has good physical-chemical properties, which are monitored in 14 wells.It is of a calcium hydrocarbonate type.The Total Hardness (TH) varies from 20 to 26 German degrees and the Total Dissolved Solids (TDS) from 650 to 900 mg/L [55,64,65].

Selection of Appropriate Groundwater Vulnerability Assessment Method
The GWV studies are continuously refined to consider the impacts of a growing population, including increased economic activities, which lead to an increase in the extracted groundwater amount and a threat to groundwater quality due to pollution [28] and overexploitation [43,66].The selection process of the GWV assessment methods and their performance evaluation should avoid contrasting findings [67].This imposes a clear understanding of the GWV notion and assessment methods and of the local factors.Moraru and Hannigan [68] think that the most appropriate definition of GWV is provided by Zaporozec, who defines it as ''an intrinsic property of a groundwater system that depends on the sensitivity of the system to human and/or natural impacts''.As such, GWV can be ''intrinsic or natural'' or ''specific or integrated" [2].Intrinsic vulnerability is a function of The groundwater extraction impact on the water level is evaluated by comparing the data obtained from the annual and monthly water level variation in three monitoring wells, located in areas of intensive groundwater extraction [61][62][63].This monitoring shows that the largest variations were recorded in the areas of intensive groundwater extraction, such as in the Mezez-Laknas and Fushe Kruje-Gjola River bridge sectors (Figure 1), where the water table varies from 2.05 m to 8.07 m and from 2.15 m to 12.2 m, respectively [63].There is no monitoring well in the recharge zones.
The natural flow from the rivers is roughly 610 L/s, while the contribution of both the net recharge and the molassic deposits is roughly 100 L/s [63].The total inflow (710 L/s) is lower than the total outflow, which varies from 1270 to 1500 L/s [25].The increase in the groundwater outflow compared to the natural inflow occurs thanks to induced filtration due to the groundwater recharging directly from the rivers and the lowering of the dynamic groundwater level.
The groundwater of the Tirana-Ishmi aquifer has good physical-chemical properties, which are monitored in 14 wells.It is of a calcium hydrocarbonate type.The Total Hardness (TH) varies from 20 to 26 German degrees and the Total Dissolved Solids (TDS) from 650 to 900 mg/L [55,64,65].

Selection of Appropriate Groundwater Vulnerability Assessment Method
The GWV studies are continuously refined to consider the impacts of a growing population, including increased economic activities, which lead to an increase in the extracted groundwater amount and a threat to groundwater quality due to pollution [28] and overexploitation [43,66].The selection process of the GWV assessment methods and their performance evaluation should avoid contrasting findings [67].This imposes a clear understanding of the GWV notion and assessment methods and of the local factors.Moraru and Hannigan [68] think that the most appropriate definition of GWV is provided by Zaporozec, who defines it as "an intrinsic property of a groundwater system that depends on the sensitivity of the system to human and/or natural impacts".As such, GWV can be "intrinsic or natural" or "specific or integrated" [2].Intrinsic vulnerability is a function of hydrogeological factors, while specific vulnerability concerns mostly the risk of the aquifer's contamination [11,69,70].Civita stresses that the specific vulnerability depends on the physical-chemical characteristics, rate of occurrence and spatial distribution of contaminants [71].Due to a lack of available data, it is impossible to assess the specific vulnerability of the Tirana-Ishmi aquifer.Thus, the intrinsic GWV is the most appropriate for the purposes of the present study.Given that there are numerous intrinsic vulnerability methods [13], we made efforts to select the most feasible of them based on the available data and information about the study area.Civita and De Maio selected more than 20 main intrinsic vulnerability methods [14][15][16]71], which are included in two main classes in terms of their use: use for any physiographic scenario or for any particular area [14,16].These classes are subdivided into three basic groups [14,16,71,72], as follows:

•
Hydrogeologic Complex and Setting assessment-HCS; Zaporozec and Vrba summarized the main worldwide research on GWV and classified vulnerability maps in terms of their purpose and content, scale and graphic representation [73].The choice of the most appropriate GWV assessment model and mapping technique should be based on a realistic evaluation of the quantity, distribution and reliability of the available data and information [13,71].This choice depends on the correlation between three main factors: the density of the monitoring points, the available information for each point and the map's scale.HCS methods and small-scale maps are used when the specific basic information is inadequate and/or scarce and scattered throughout the area under study.Parametric System assessment and medium-scale maps are used when there is a moderate monitoring point density, with a satisfying spatial distribution and information, while large-scale maps are prepared when the available information per unit area is sufficient [16,71].AR and numerical model assessment take into consideration only a part of the features of the area to be assessed.The AR methods proposed by Zampetti [74] and Fried [75] take into consideration only the unsaturated zone and water table depth, to which Marcolongo and Pretto [76] added the net recharge and the effective soil moisture.As the data reliability in an area can widely vary with the mean elevation a.s.l., simple methods such as HCS and MS should be used for mapping the vulnerability of hilly and mountainous areas.By contrast, the parametric RS and PCSM methods are mostly recommended for flat and foothill areas, when the amount and reliability of data, measurements, tests and analyses can be considered to be sufficient for the required mapping scale [16,71].The parametric methods are an extension of the HCS methods, as they use a rating or a numerical score for each hydrogeological parameter included and judged necessary and adequate for assessing GWV.In addition to the rating score, the PCSM methods assign an importance weight to each considered parameter, in order to reflect the relationships between the different parameters, as well as their relevance in the vulnerability assessment.PCSM methods can be easily implemented, modified, updated and used for large areas [69].
Based on the above and the available data and information on the Tirana-Ishmi aquifer, the PCSM methods' results are the most suitable for the assessment of this aquifer.Civita and De Maio [16,71] distinguish four main PCSM methods; DRASTIC (Aller and Thornhill, 1987) [77], Trojan and Perry (1988) [78], ISIS (Civita and De Regibus, 1995) [79] and SINTACS (Civita, 1990 [80] and 1994 [14]; and Civita and De Maio, 1997 [15]).The Trojan and Perry method is used mostly for areas of large extensiveness [78].The ISIS method is applied in both porous and karst aquifer systems [79].The Trojan and Perry and ISIS methods do not take into account the aquifer's hydraulic conductivity [13,16].Amongst the four main PCSM methods, the most used worldwide are DRASTIC and SINTACS, which take into consideration more parameters than the other ones.Although DRASTIC remains the most popular [69,81], the ratings and weights of the factors taken into consideration by this method need to be adapted to the local features [69,82,83].Although both the DRASTIC and SINTACS models use the same parameters, the rating and weighting of the SINTACS procedure is more flexible to specific features of a study area [8,69,84,85].
Gogu and Dassargues prefer the methods resulting in a more contrasted sensitivity for a specific area, because their results can be used and interpreted better [86].This preference can be satisfied using SINTACS, which allows one to assign different weighting rates to different zones of the same aquifer in terms of the local features.Similarly to other PCSM methods, the results of SINTACS can be influenced by subjectivity in the weighting rate assignment [69].However, this subjectivity can be reduced through the statistical analysis and validation of the used method [87,88].SINTACS was created, elaborated and continuously improved by Italian researchers [14][15][16][17]71,80].It is a development of the DRASTIC model adapted to the geological, hydrogeological, soil and climatic conditions that are typical to Italy and therefore to other Mediterranean countries, including Albania, whose natural features are more or less comparable to those of Italy.Based on the characteristics of the study area, the available data and information and the appropriate scale of the vulnerability map and its purpose, as well as on the characteristics of SINTACS model [6,89], SINTACS stands out as the most appropriate method for assessing the vulnerability of the Tirana-Ishmi aquifer.Furthermore, as the study area is a flat alluvial plain [21,22], the reliability of the vulnerability assessment using SINTACS should be high [14][15][16]90].The map's scale, 1:50,000, is preferred for the purposes of the present study.The use of a GIS for mapping the aquifer's vulnerability allows for the continuous update and improvement of the vulnerability assessment whenever the data and information are updated and/or completed.
Each parameter is assigned a range of values divided into discrete intervals from 1 to 10.Each interval is given a rate which reflects the relative vulnerability degree.Each parameter is multiplied by a weighting coefficient, which expresses the relationships amongst the parameters and their relevance to the GWV assessment.The rate given to each interval is multiplied by the weighting string of each parameter, and the products are summed, thus giving the vulnerability index (VI), as provided by Formula (1).
where P i -parameter's value and W i -relative weight of the selected weighting strings.
Based on the vulnerability index, the study area is split into polygons.The interval of the vulnerability index values for each polygon ranges from 26 to 260 [15].The higher the vulnerability index, the more sensitive to pollution is the groundwater system.The SINTACS model uses five weighting strings which correspond to normal impact, high impact, with drainage, karst and fissured rocks [15].Based on the geological and hydrogeological settings of the study area, and on its socio-economic aspects, the most suitable weighting string corresponds to normal impact.Table 1 provides the theoretical weighs for normal impact, in accordance with the provisions of the SINTACS method [92].
The vulnerability degree is divided into classes in terms of the vulnerability index, as provided in Table 2 [92].

Tools for Aquifer's Vulnerability Mapping
Each SINTACS parameter map and the GWV map were prepared through ArcGIS software [93], while the aquifer 3D model was created based on the interpolation of 35 water wells' data, using the Groundwater Modeling System (GMS), version 10.7 [59].

Validation
The validation of both the selected vulnerability model and the prepared vulnerability map serves to avoid any eventual subjectivity that may lead to an inaccurate assessment [86].The validation was performed through a sensitivity analysis (SA), which expresses the effect of change in one factor on another one.SA influences the model formulation, calibration and verification [19], and helps to analyze the interaction between parameters, their appropriate range and spatial variability, and to identify the parameters having the highest impact on the model response and outputs, which in turn influence the model outcomes [94].

Sensitivity Analysis
The selected SA approach includes the map removal SA and the single-parameter SA [87,88,[95][96][97].The map removal SA expresses the effect of each parameter's map on the vulnerability index [87].Formula (2) provides the SA in terms of the variation index Si.

Si
where V i is the total vulnerability index for N = 7 layers, while V pi is the vulnerability index corresponding to n = 6 layers after the removal of the parameter p map.The single-parameter SA helps to compare the effective weights for each parameter to their theoretical weights [88] provided by the SINTACS method.Formula (3) calculates the effective weight for each polygon of the vulnerability map.
where P ri and P wi are, respectively, the rates and the weight of the layer P of the polygon i, while VI is the vulnerability index of the polygon i.

Validation Using the Nitrates Content in Groundwater
The accuracy of the selected vulnerability assessment method was validated using nitrates content (NO 3 in mg/L), which is an anthropogenic pollution indicator.The NO 3 concentration in groundwater is largely used and recommended worldwide for validating vulnerability assessment methods [84,[97][98][99].The nitrates content in the groundwater comes from different contaminants [100], including fertilizers and the municipal wastewater that flows into the local rivers and constitutes the main source of groundwater pollution in the study area [101][102][103].The data on NO 3 were taken from the State of Environment Report (SoER), which is annually published by the National Environment Agency (NEA) [104].The monitoring period extents from 2015 to 2022.

Parameters Maps
The seven SINTACS parameters [92] were discussed and mapped based on the available data and information.
Parameter "S": the water table depth.The water table depth is in inverse proportion to the natural vulnerability.The assigned SINTACS rates to this parameter decrease with the increase in the water table depth, following a hyperbolic curve.The assigned range varies from 10 to 1.The rate 10 is assigned to the recharge zones, and when the water table is on the ground surface.SINTACS rates should be corrected for confined aquifers, especially for artesian ones [92].The data on water levels in the wells are collected from hydrogeological monitoring projects carried out by the AGS during the 2015-2022 period [104] and from the hydrogeological map [21], as well as from several published studies [49,50,[53][54][55].
In the Tirana City area and the terraces of the Tirana, Terkuza, Zeza and Droja Rivers, the aquifer is unconfined (or phreatic).The water wells document that the cover layer thickness in these sectors varies from 0.3 m to 5m, with the exception of the river beds, where the cover layer is nonexistent.The water table depth varies from zero to 5.0 m from the ground surface and is subject to seasonal precipitation [49,53,54].The SINTACS rate assigned to parameter "S" in these sectors is 9 [92] (Table 3).In the Laknas-Berxulle area that lies north of Tirana City, the aquifer was confined until about 20 years ago [49,50,53], but nowadays it has changed to unconfined due to overexploitation [23,104].As the water table depth in this area varies from 5 m to 10 m, the assigned SINTACS rate is 7 [92].In the southern edge of Tirana and in the sector extending from the central part of the study area towards its NW, the aquifer is confined.According to Civita, when an aquifer is confined, the water table depth is considered equal to the cover layer's thickness [92].In the SINTACS range assessment, the assigned rates corresponding to cover layer thicknesses of 10-20 m, 20-30 m and greater than 30 m are 4, 2.5 and 2-1, respectively (Table 3).For water table thicknesses greater than 30 m, the assigned rate is 1.Table 3 provides the depth intervals and the evaluation of the "S" parameter, based on the water level depth analysis.
The water table depth parameter map (Figure 5) shows that roughly up to 65% of the study area's surface is assigned the lowest rates (1, 2 and 4) and up to 25%, rate 7, whereas the remaining 10% earns the highest rate (9).The lowest rate belongs to the confined parts of the aquifer, where the piezometric surface is close to the ground surface, whereas the highest rate is associated with the recharge zones located in the rivers' terraces.
The water table depth parameter map (Figure 5) shows that roughly up to 65 study area's surface is assigned the lowest rates (1, 2 and 4) and up to 25%, rate 7, the remaining 10% earns the highest rate (9).The lowest rate belongs to the confin of the aquifer, where the piezometric surface is close to the ground surface, whe highest rate is associated with the recharge zones located in the rivers' terraces.Parameter "I": the net recharge.Precipitation is the main aquifer recharge Basic information on the net recharge (or effective infiltration) is obtained from a continuous data collected in the period 1961-2020 from the former Institute of Hy teorology (HIDMET) [41] and nowadays from the Institute of Geosciences (IGEO) 10-15% of the precipitation amount infiltrates through the cover layer, while the ing goes to evapotranspiration and surface runoff [41].The rainfalls' distribution i uniform for the whole Albanian Western Lowland where the study area is includ The net recharge is calculated according to Turc [105] for an average precipitation mm/year and an average temperature of 15 °C.The net recharge is calculated as th uct of the infiltration with the potential infiltration coefficient χ.Table 4 provides tential infiltration coefficient of a cover layer consisting of medium to fine alluvial [92].
In the areas where the cover layer thickness varies from zero to 5.0 m, the p infiltration coefficient is χ = 0.3, and the calculated effective infiltration is 213 m The assigned SINTACS rate is 9.In the aquifer sectors where the cover layer th varies from 5 m to 10 m and the lithological composition includes clays, sand an with a potential infiltration coefficient χ = 0.2, a rate 6 is assigned (Table 4).In the with a cover layer of 10 m to 20 m thickness, the assigned SINTACS rate is 5.By c in the confined part of the aquifer in the northwest of the study area, due to the h content and large thickness of the cover layer, the impact of the net recharge is pr nonexistent, and therefore the rate zero has been assigned to the parameter "I".Parameter "I": the net recharge.Precipitation is the main aquifer recharge source.Basic information on the net recharge (or effective infiltration) is obtained from a series of continuous data collected in the period 1961-2020 from the former Institute of Hydrometeorology (HIDMET) [41] and nowadays from the Institute of Geosciences (IGEO).Nearly 10-15% of the precipitation amount infiltrates through the cover layer, while the remaining goes to evapotranspiration and surface runoff [41].The rainfalls' distribution is almost uniform for the whole Albanian Western Lowland where the study area is included [41].The net recharge is calculated according to Turc [105] for an average precipitation of 1250 mm/year and an average temperature of 15 • C. The net recharge is calculated as the product of the infiltration with the potential infiltration coefficient χ.Table 4 provides the potential infiltration coefficient of a cover layer consisting of medium to fine alluvial deposits [92].In the areas where the cover layer thickness varies from zero to 5.0 m, the potential infiltration coefficient is χ = 0.3, and the calculated effective infiltration is 213 mm/year.The assigned SINTACS rate is 9.In the aquifer sectors where the cover layer thickness varies from 5 m to 10 m and the lithological composition includes clays, sand and gravel with a potential infiltration coefficient χ = 0.2, a rate 6 is assigned (Table 4).In the sectors with a cover layer of 10 m to 20 m thickness, the assigned SINTACS rate is 5.By contrast, in the confined part of the aquifer in the northwest of the study area, due to the high clay content and large thickness of the cover layer, the impact of the net recharge is practically nonexistent, and therefore the rate zero has been assigned to the parameter "I".
Figure 6 shows the net recharge parameter map, which indicates that the assigned SINTACS rate for this parameter is zero in roughly 50% of the study area, where the aquifer is confined and the cover layer is thick and composed of fine alluvial deposits.In 40% of the study area, this parameter has the medium rate, 5-6, whereas the remaining 10% has high rates, assigned to the rivers' terraces, where the cover layer is nonexistent or is thinner than 5 m and non-uniformly spread.Figure 6 shows the net recharge parameter map, which indicates that the assig SINTACS rate for this parameter is zero in roughly 50% of the study area, where the uifer is confined and the cover layer is thick and composed of fine alluvial deposit 40% of the study area, this parameter has the medium rate, 5-6, whereas the remai 10% has high rates, assigned to the rivers' terraces, where the cover layer is nonexiste is thinner than 5 m and non-uniformly spread.Parameter "N": unsaturated media.Both the soil media and unsaturated media tect the aquifer from pollution, thanks to their physical and chemical properties.The ture of the unsaturated media determines the travel time of a contaminant [92].The li ogy of the unsaturated zone observed in situ or in the drilled wells [23,24,49,50] play important role in the creation of a conceptual model of the hydrogeological system of zone.In the 3D model created for this study's purposes (Figure 4), the lithological are correlated.The reliability of this model was verified based on the lithotypes prov in the geotechnical maps [106][107][108][109][110][111].Based on the lithology of each polygon, the valu parameter "N" was defined using the weighted average calculation provided in Form (4).Parameter "N": unsaturated media.Both the soil media and unsaturated media protect the aquifer from pollution, thanks to their physical and chemical properties.The texture of the unsaturated media determines the travel time of a contaminant [92].The lithology of the unsaturated zone observed in situ or in the drilled wells [23,24,49,50] plays an important role in the creation of a conceptual model of the hydrogeological system of this zone.In the 3D model created for this study's purposes (Figure 4), the lithological data are correlated.The reliability of this model was verified based on the lithotypes provided in the geotechnical maps [106][107][108][109][110][111].Based on the lithology of each polygon, the value of parameter "N" was defined using the weighted average calculation provided in Formula (4).

N =
where P i is the SINTACS rate of each lithotype (Table 5), h i is the thickness of each lithotype and h is the unsaturated zone's thickness.

SINTACS Range
Coarse alluvial deposits 8-9 Fine to medium alluvial deposits 6-8 Clay, silt 1-3 Table 5 provides the lithotypes of the unsaturated zone and the respective SINTACS range corresponding to alluvial deposits [92], while Figure 7 shows the map of the unsaturated zone parameter.
where Pi is the SINTACS rate of each lithotype (Table 5), hi is the thickness of each l type and h is the unsaturated zone's thickness.
Table 5 provides the lithotypes of the unsaturated zone and the respective SINT range corresponding to alluvial deposits [92], while Figure 7 shows the map of the saturated zone parameter.

SINTACS
Range Coarse alluvial deposits 8-9 Fine to medium alluvial deposits 6-8 Clay, silt 1-3 According to Figure 7, the unsaturated zone parameter rates indicate a larger sp variability compared to the other parameters.The highest rates are associated with river terraces (rate 9) and Tirana City area (rate 7) that constitute roughly 10% and 20 the study area's surface, respectively.Roughly up to 30% of the study area, belongin the central and northeastern parts of this area, is assigned an almost medium rate (ra The lowest rates (1, 2 and 3) belong to about 40% of the study area, mainly in its we and northwestern parts, where the unsaturated zone is composed of thick clays and Parameter "T": soil media.The soil media are the first barrier to groundwater p tion.The soil classification is based on the World Reference Base for Soil Resources (W [112,113].On the basis of their physical-morphological characteristics, the soil hori are distinguished from their composition, the grain diameter size and the texture structure [113][114][115][116][117]. A large number of physical-chemical analyses have been effectu by the Albanian Institute of Soil Science in the framework of different scientific pro whose findings were synthesized in the preparation of a geo-referenced soil databas Albania [118,119], followed by soil mapping at 1:250,000 scale for the whole country According to Figure 7, the unsaturated zone parameter rates indicate a larger spatial variability compared to the other parameters.The highest rates are associated with the river terraces (rate 9) and Tirana City area (rate 7) that constitute roughly 10% and 20% of the study area's surface, respectively.Roughly up to 30% of the study area, belonging to the central and northeastern parts of this area, is assigned an almost medium rate (rate 6).The lowest rates (1, 2 and 3) belong to about 40% of the study area, mainly in its western and northwestern parts, where the unsaturated zone is composed of thick clays and silts.
Parameter "T": soil media.The soil media are the first barrier to groundwater pollution.The soil classification is based on the World Reference Base for Soil Resources (WRB) [112,113].On the basis of their physical-morphological characteristics, the soil horizons are distinguished from their composition, the grain diameter size and the texture and structure [113][114][115][116][117]. A large number of physical-chemical analyses have been effectuated by the Albanian Institute of Soil Science in the framework of different scientific projects, whose findings were synthesized in the preparation of a geo-referenced soil database for Albania [118,119], followed by soil mapping at 1:250,000 scale for the whole country and at 1:50,000 scale for the Western Lowland of Albania [120,121].According to the WRB classification, the soil types within the study area are of the Fluvisol, Cambisol, Arenosol and Luvisol types [121] (Figure 8).Table 6 provides the study area's soil types and the respective SINTACS rates [92].Figure 8 shows the soil media parameter map, which indicates that roughly up to 5% of the aquifer has a high rate ( 8), 88% has rate 5 and 4.5 and 7% has a low rate (2.5).
Hydrology 2024, 11, x FOR PEER REVIEW 15 at 1:50,000 scale for the Western Lowland of Albania [120,121].According to the WRB sification, the soil types within the study area are of the Fluvisol, Cambisol, Arenosol Luvisol types [121] (Figure 8).Table 6 provides the study area's soil types and the res tive SINTACS rates [92].Figure 8 shows the soil media parameter map, which indicates that roughly up t of the aquifer has a high rate ( 8), 88% has rate 5 and 4.5 and 7% has a low rate (2.5).Parameter "A": aquifer media.When a pollutant penetrates through both the media and the unsaturated zone and reaches the aquifer, its impact can be mitig through the natural processes of dilution, dispersion and absorption, as well as by chemical processes that may occur in the geological formation.
The influence of the aquifer media characteristics on GWV includes the followin  Parameter "A": aquifer media.When a pollutant penetrates through both the soil media and the unsaturated zone and reaches the aquifer, its impact can be mitigated through the natural processes of dilution, dispersion and absorption, as well as by the chemical processes that may occur in the geological formation.
The influence of the aquifer media characteristics on GWV includes the following: • The recharge zones, which are located in the riverbeds in the eastern edge of the study area; • The unconfined areas of the aquifer; • The confined areas of the aquifer; • The heterogeneity of the aquifer's particles size; • The groundwater flow direction.
Based on the available data and information on the above aquifer media's features [21,23,24,54,57,58], Table 7 provides the assigned SINTACS rate, in terms of the lithotypes of the saturated zone [92].Figure 9 shows the aquifer media parameter map, which indicates that the highest rate (9 and 8) is assigned to roughly 50% of the aquifer, located in the southern, central and eastern parts of the study This high rate is due to the presence of the recharge zones, the unconfined aquifer sections and the coarse grain-size alluvium deposits.The latter are found mainly in the southern and eastern part of the study area, close to the foothills.In the western and northwestern part of the study area, the assigned SINTACS rate for this parameter is 6 and 7, corresponding to 20% and 30% of the area, respectively.
Hydrology 2024, 11, x FOR PEER REVIEW Based on the available data and information on the above aquifer media's [21,23,24,54,57,58], Table 7 provides the assigned SINTACS rate, in terms of the lit of the saturated zone [92].Figure 9 shows the aquifer media parameter map, which indicates that the rate (9 and 8) is assigned to roughly 50% of the aquifer, located in the southern and eastern parts of the study area.This high rate is due to the presence of the r zones, the unconfined aquifer sections and the coarse grain-size alluvium depo latter are found mainly in the southern and eastern part of the study area, clos foothills.In the western and northwestern part of the study area, the assigned S rate for this parameter is 6 and 7, corresponding to 20% and 30% of the area, resp Parameter "C": hydraulic conductivity.The hydraulic conductivity controls of the groundwater flow through the porous media.The higher the hydraulic co ity, the higher is the aquifer's vulnerability.The assigned SINTACS rate to this pa for each polygon (Table 8) is based on the experimental data on the hydraulic cond of the study area [49,50,53,54].Parameter "C": hydraulic conductivity.The hydraulic conductivity controls the rate of the groundwater flow through the porous media.The higher the hydraulic conductivity, the higher is the aquifer's vulnerability.The assigned SINTACS rate to this parameter for each polygon (Table 8) is based on the experimental data on the hydraulic conductivity of the study area [49,50,53,54].Figure 10 shows the hydraulic conductivity parameter map, which indicates that the highest hydraulic conductivity rates (8-9) belong to roughly 50% of the aquifer, while rate 7 is assigned to the northwestern part of the study area, where the grain size of the sediments composing the aquifer is smaller Figure 10 shows the hydraulic conductivity parameter map, which indicates highest hydraulic conductivity rates (8-9) belong to roughly 50% of the aquifer, w 7 is assigned to the northwestern part of the study area, where the grain size of t ments composing the aquifer is smaller Parameter "S": topographic slope.The topography controls whether the po will run away or will be infiltrated throughout the cover layer and reach and a groundwater.The terrain slope map of the study area is based on the topographic slope of each elementary polygon.According to the digital terrain model map, wh tical accuracy is ±10 cm for urban areas, ±20 cm for rural areas and ±40 cm for moun areas [36], the study area is located mainly on flat terrain (slope of 0-2%), with th tion of a few small areas with a very low slope of 2-4%.By classifying the who area at a slope of 0-4%, the assigned SINTACS rate to the parameter S is 9 [92].F shows the terrain and the topographic slope parameter maps.Parameter "S": topographic slope.The topography controls whether the pollutants will run away or will be infiltrated throughout the cover layer and reach and affect the groundwater.The terrain slope map of the study area is based on the topographic surface slope of each elementary polygon.According to the digital terrain model map, whose vertical accuracy is ±10 cm for urban areas, ±20 cm for rural areas and ±40 cm for mountainous areas [36], the study area is located mainly on flat terrain (slope of 0-2%), with the exception of a few small areas with a very low slope of 2-4%.By classifying the whole study area at a slope of 0-4%, the assigned SINTACS rate to the parameter S is 9 [92].Figure 11 shows the terrain and the topographic slope parameter maps.

Statistical Analysis of the SINTACS Parameters
To interpret their reliability, the obtained results were submitted to a statistical analysis [18], which helps also to compare the relevance of each of the seven parameters in the vulnerability index [122].The mean value, the standard deviation and the variation coefficient were analyzed for each parameter [123,124].The standard deviation measures the amount of variation or dispersion in a set of data values relative to its mean, while the coefficient of variation, also called the relative standard deviation, is the ratio of the standard deviation to the mean value.The lower the standard deviation and the coefficient of variation, the closer to the mean are the parameters' values.Table 9 provides statistical data on the seven parameters used in the SINTACS method for vulnerability mapping.9 indicates that the unsaturated zone has the highest impact on the vulnerability index, as its mean has the highest value (24.37), followed by hydraulic conductivity (23.36) and the aquifer media (23.04).The soil media and the water table depth show al-

Statistical Analysis of the SINTACS Parameters
To interpret their reliability, the obtained results were submitted to a statistical analysis [18], which helps also to compare the relevance of each of the seven parameters in the vulnerability index [122].The mean value, the standard deviation and the variation coefficient were analyzed for each parameter [123,124].The standard deviation measures the amount of variation or dispersion in a set of data values relative to its mean, while the coefficient of variation, also called the relative standard deviation, is the ratio of the standard deviation to the mean value.The lower the standard deviation and the coefficient of variation, the closer to the mean are the parameters' values.Table 9 provides statistical data on the seven parameters used in the SINTACS method for vulnerability mapping.Table 9 indicates that the unsaturated zone has the highest impact on the vulnerability index, as its mean has the highest value (24.37), followed by hydraulic conductivity (23.36) and the aquifer media (23.04).The soil media and the water table depth show almost the same mean value, 21.55 and 20.97, respectively.The net recharge has the minimum mean value (14.28).The very high value (94.42) of the variation coefficient related to net recharge indicates that the vulnerability index is highly influenced by this parameter, even when it varies slightly.Although the soil media parameter and the water table depth have more or less the same mean value, their variation coefficients are different-respectively, 31.93 and 70.72-indicating that the water table depth has a higher contribution to the variation in the vulnerability index than the soil media.Lastly, because of the zero value of its variation coefficient, it is not expected that the topographic slope will influence the vulnerability index variation.

Vulnerability Map According to the Theoretical Weights of the SINTACS Parameters
Table 10 provides the vulnerability classes for the Tirana-Ishmi aquifer, and the related surface areas for each class, in compliance with the SINTACS assessment [92], using GIS mapping [14][15][16][17].It indicates that, based on the vulnerability index, there are four vulnerability classes within the Tirana-Ishmi aquifer: low, moderate, high and extremely high.The GWV for normal impact is mapped using the GIS (Figure 12).As it has not yet been submitted to a sensitivity analysis, the prepared map is called a preliminary vulnerability map.mum mean value (14.28).The very high value (94.42) of the variation coefficient related net recharge indicates that the vulnerability index is highly influenced by this parame even when it varies slightly.Although the soil media parameter and the water table de have more or less the same mean value, their variation coefficients are different-resp tively, 31.93 and 70.72-indicating that the water table depth has a higher contribution the variation in the vulnerability index than the soil media.Lastly, because of the z value of its variation coefficient, it is not expected that the topographic slope will influe the vulnerability index variation.

Vulnerability Map According to the Theoretical Weights of the SINTACS Parameters
Table 10 provides the vulnerability classes for the Tirana-Ishmi aquifer, and the lated surface areas for each class, in compliance with the SINTACS assessment [92], us GIS mapping [14][15][16][17].It indicates that, based on the vulnerability index, there are fo vulnerability classes within the Tirana-Ishmi aquifer: low, moderate, high and extrem high.The GWV for normal impact is mapped using the GIS (Figure 12).As it has not been submitted to a sensitivity analysis, the prepared map is called a preliminary vuln ability map.The preliminary vulnerability map shows that the extremely high vulnerability class is associated with the recharge zones of the Tirana, Terkuza, Gjola and Zeza riverbeds, within the eastern part of the study area.The high vulnerability class characterizes the southern part of the study area where the aquifer is mostly unconfined, as well as some sections in the study area eastern part, where the cover layer consists mostly of medium grain-size alluvial deposits.The moderate vulnerability class characterizes mainly some confined aquifer parcels that border from the east the aquifer section of low vulnerability.The latter includes the confined aquifer sections, which are overlaid by a thick clayey cover layer.

Validation of Preliminary Vulnerability Map Using Sensitivity Analysis
The validation of the preliminary vulnerability map was performed using the SA that includes the map removal sensitivity analysis and the single-parameter SA.In addition, the NO 3 content measured in situ is taken into consideration as an indicator of anthropogenic pollution [84,[97][98][99][100].The SA approach used in this study includes the map removal SA and the single-parameter SA.
Map removal sensitivity analysis.The aim of the map removal SA is to evaluate the relevance of each of the seven SINTACS parameters in the vulnerability assessment.Lodwick et al. [87] suggest to analyze the sensitivity by removing one parameter map, to better understand the relevance of each parameter (Table 11).The statistical parameters of the variation index in Table 11 are calculated based on Equation (2), by each time removing one parameter map.Table 11 indicates that the most significant parameter in the variation index (SI) is the net recharge (mean value = 1.13), which has the highest standard deviation value.That is probably related to the high value (4) of its theoretical weight.The variation index is also sensitive to the removal of the unsaturated zone, followed by the depth to groundwater.The mean values of these parameters are respectively 1.02 and 0.87.This may result from their high theoretical weights, which, according to the SINTACS model, have a value of 5 for both these parameters.The relatively high impact of the unsaturated zone is associated with the cover layer's influence on aquifer vulnerability.Similarly, the water table depth has a relatively high impact, but slightly lower than that of the unsaturated zone, because the aquifer is mainly confined, whereas the aquifer media have less effect on the variation index, which shows the lowest mean value (0.53) after the removal of this parameter.This is probably due to the increase in the sandy fraction towards the northwest of the aquifer.Likewise, the variation index shows a low mean value (0.54) after the removal of the topographic slope map.This low value is linked to the low (2) theoretical weight.Since their mean values are 0.67 and 0.60, respectively, hydraulic conductivity and soil media have a moderate contribution to the variation index.
Single-parameter sensitivity analysis.This analysis compares the parameters' effective weights with their respective SINTACS model's theoretical weights (Table 12). Figure 13 shows the mean effective weights and the mean theoretical weights of the SINTACS parameters.
Single-parameter sensitivity analysis.This analysis compares the parameters' effective weights with their respective SINTACS model's theoretical weights (Table 12). Figure 13 shows the mean effective weights and the mean theoretical weights of the SINTACS parameters.13 indicate that the effective weight of the water table depth (S) is 36% lower than its theoretical value, which seems to be overestimated.This difference may be due to the fact that in confined aquifers, the water table depth is not a significant parameter.Therefore, when the aquifer is confined, the theoretical value of the water table depth is lowered.The net recharge (I) has an effective weight twice lower than its theoretical one.Table 11 indicates that this parameter has the greatest impact in the vulnerability index.This can be due to the overestimation of its theoretical weight.The net recharge is low or negligible in roughly 70% of the study area.This parameter is linked to the thickness and lithological composition of the cover layer.The mean effective weights of the unsaturated zone (N) and the soil media (T) are, respectively, nearly 14% and 7% lower than their theoretical values.Therefore, their effective and theoretical weights are closer to each other compared to the other SINTACS parameters.The unsaturated zone and the soil media constitute, respectively, the second and first barrier to groundwater pollution Both the mean effective weights of the aquifer media (A) and hydraulic conductivity (C) are 50% higher than their theoretical ones.Because of the flat terrain, the topographic slope (S) has a mean effective weight nearly 1.8 times greater than its SINTACS theoretical value.
Based on the above, the aquifer media, hydraulic conductivity and topographic slope seem to be underestimated in the selected vulnerability assessment model; whereas, the   13 indicate that the effective weight of the water table depth (S) is 36% lower than its theoretical value, which seems to be overestimated.This difference may be due to the fact that in confined aquifers, the water table depth is not a significant parameter.Therefore, when the aquifer is confined, the theoretical value of the water table depth is lowered.The net recharge (I) has an effective weight twice lower than its theoretical one.Table 11 indicates that this parameter has the greatest impact in the vulnerability index.This can be due to the overestimation of its theoretical weight.The net recharge is low or negligible in roughly 70% of the study area.This parameter is linked to the thickness and lithological composition of the cover layer.The mean effective weights of the unsaturated zone (N) and the soil media (T) are, respectively, nearly 14% and 7% lower than their theoretical values.Therefore, their effective and theoretical weights are closer to each other compared to the other SINTACS parameters.The unsaturated zone and the soil media constitute, respectively, the second and first barrier to groundwater pollution.Both the mean effective weights of the aquifer media (A) and hydraulic conductivity (C) are 50% higher than their theoretical ones.Because of the flat terrain, the topographic slope (S) has a mean effective weight nearly 1.8 times greater than its SINTACS theoretical value.
Based on the above, the aquifer media, hydraulic conductivity and topographic slope seem to be underestimated in the selected vulnerability assessment model; whereas, the water table depth and the net recharge are overestimated.Regarding the unsaturated zone and soil media, the model weights are close to the effective weights.Therefore, for a realistic parameters assessment, accurate, detailed and representative information is needed for all these parameters.

Vulnerability Map Using Effective Parameters' Weights
The last step of GWV mapping is the preparation of the vulnerability map using effective weights, which hereinafter is called the final vulnerability map or, simply, the vulnerability map (Figure 14).
Hydrology 2024, 11, x FOR PEER REVIEW 22 water table depth and the net recharge are overestimated.Regarding the unsaturated and soil media, the model weights are close to the effective weights.Therefore, for a istic parameters assessment, accurate, detailed and representative information is ne for all these parameters.

Vulnerability Map Using Effective Parameters' Weights
The last step of GWV mapping is the preparation of the vulnerability map usin fective weights, which hereinafter is called the final vulnerability map or, simply, the nerability map (Figure 14).The (final) GWV map includes five vulnerability classes: "extremely high", " high", "high", "moderate" and "low".The comparison of this map (Figure 14) to the oretical weights' map (or preliminary map) (Figure 12), indicates that the areas of an tremely high" vulnerability class are the same in both maps and represent only th charge zones.The "very high" vulnerability class is not present in the preliminary m This class appears in the final map, where it is located in the areas of the "high vuln bility" class of the preliminary map. Figure 14 shows that the "very high" vulnerab class is also encountered in small areas located mainly westward of the extremely vulnerability class areas.In the areas of a "very high" vulnerability class, the aquif mostly unconfined or has been transformed into unconfined due to groundwater ove ploitation.This occurs mainly in the southern part of the study area, which is highl banized.Compared to the preliminary map, in the final map, there is a noticeable decr in the surface areas with "high", "moderate" and "low" vulnerability classes, as well displacement of these areas northwestward and westward.Table 13 indicates the vu ability classes of the Tirana-Ishmi aquifer based on the effective weights, as well as spatial extent.The (final) GWV map includes five vulnerability classes: "extremely high", "very high", "high", "moderate" and "low".The comparison of this map (Figure 14) to the theoretical weights' map (or preliminary map) (Figure 12), indicates that the areas of an "extremely high" vulnerability class are the same in both maps and represent only the recharge zones.The "very high" vulnerability class is not present in the preliminary map.This class appears in the final map, where it is located in the areas of the "high vulnerability" class of the preliminary map. Figure 14 shows that the "very high" vulnerability class is also encountered in small areas located mainly westward of the extremely high vulnerability class areas.In the areas of a "very high" vulnerability class, the aquifer is mostly unconfined or has been transformed into unconfined due to groundwater overexploitation.This occurs mainly in the southern part of the study area, which is highly urbanized.Compared to the preliminary map, in the final map, there is a noticeable decrease in the surface areas with "high", "moderate" and "low" vulnerability classes, as well as a displacement of these areas northwestward and westward.Table 13 indicates the vulnerability classes of the Tirana-Ishmi aquifer based on the effective weights, as well as their spatial extent.Apart the validation performed using the sensitivity analysis, the reliability of the selected vulnerability assessment model and of the prepared final vulnerability map were also validated using the nitrates content, which is an indicator of the groundwater's anthropogenic pollution [84,[97][98][99][100].However, it should be noted that the use of nitrates content has some limitations, due to the insufficient number of monitoring wells in the study area and their spatial distribution.As shown in the vulnerability map (Figure 14), the 14 monitoring wells are not uniformly distributed in the areas associated with different vulnerability classes.In particular, there are no monitoring wells in the areas with "extremely high" and "low" vulnerability classes.In the areas with "very high", "high" and "moderate" vulnerability, the number of monitoring wells is six, five and three, respectively, which is not sufficient.
The accuracy of the preliminary vulnerability map is evaluated using the nitrates content as a validation indicator (Figure 15).Although the SINTACS parameters were not yet submitted to sensitivity analysis, Figure 15 indicates a satisfying correlation between the vulnerability classes and the variation trend of the nitrates content in the groundwater.Therefore, it can be said that the use of the nitrates content as a validation technique supports the selection of SINTACS as a GWV model and the necessity of the performed statistical analysis.In Figure 15, the average NO 3 content in the monitoring wells for the period 2015-2022 is, respectively, 9.45 mg/L in the "low" vulnerability class, 11.65 mg/L in the "moderate" vulnerability class and 18.3 mg/L in the "high" vulnerability class.The highest NO 3 content is 36 mg/L, which is lower than the maximum value allowed by national standards (50 mg/L).It should be noted that the Albanian drinking water quality standards [125] are identical to the EU Directive 98/83/EC standards [126], which have been updated by EU Directive 2020/2184 [127].

Validation Using Nitrates Content in Groundwater
Apart the validation performed using the sensitivity analysis, the reliability of the selected vulnerability assessment model and of the prepared final vulnerability map were also validated using the nitrates content, which is an indicator of the groundwater's anthropogenic pollution [84,[97][98][99][100].However, it should be noted that the use of nitrates content has some limitations, due to the insufficient number of monitoring wells in the study area and their spatial distribution.As shown in the vulnerability map (Figure 14), the 14 monitoring wells are not uniformly distributed in the areas associated with different vulnerability classes.In particular, there are no monitoring wells in the areas with "extremely high" and "low" vulnerability classes.In the areas with "very high", "high" and "moderate" vulnerability, the number of monitoring wells is six, five and three, respectively, which is not sufficient.
The accuracy of the preliminary vulnerability map is evaluated using the nitrates content as a validation indicator (Figure 15).Although the SINTACS parameters were not yet submitted to sensitivity analysis, Figure 15 indicates a satisfying correlation between the vulnerability classes and the variation trend of the nitrates content in the groundwater.Therefore, it can be said that the use of the nitrates content as a validation technique supports the selection of SINTACS as a GWV model and the necessity of the performed statistical analysis.In Figure 15, the average NO3 content in the monitoring wells for the period 2015-2022 is, respectively, 9.45 mg/L in the "low" vulnerability class, 11.65 mg/L in the "moderate" vulnerability class and 18.3 mg/L in the "high" vulnerability class.The highest NO3 content is 36 mg/L, which is lower than the maximum value allowed by national standards (50 mg/L).It should be noted that the Albanian drinking water quality standards [125] are identical to the EU Directive 98/83/EC standards [126], which have been updated by EU Directive 2020/2184 [127].Figure 16 shows that the higher the NO3 content in the groundwater, the higher is the aquifer vulnerability and vice versa.The overall low level of NO3 content in the whole aquifer could be due to the very low use of fertilizers in agriculture during the last 30 years.There are no studies on nitrates' travel time throughout the cover layer in the study area.The travel time of the nitrates depends on the water table depth, the hydraulic properties of the cover layer, the soil type, etc. [101].Bancheri et al. calculated that the average travel time of the nitrates through a clayey-rich cover layer and the overlaid soil in a case study in Italy was roughly 37 years [101].The higher the clayey content in the cover layer and soil, the longer is the travel time of the nitrates from the surface to the groundwater.Regarding the Tirana-Ishmi aquifer, the highest nitrates content is observed in the Tirana City area and its neighboring area, where the aquifer is phreatic, the water table is almost shallow, the cover layer is composed of medium grain-size alluvium deposits, and both the levels of the groundwater extraction rate and the amount of wastewater discharged into the local rivers' waters are high.The sources of nitrates in groundwater in this part of the study area are mostly sewage and industrial activities [29,100,103,128].The lowest level of NO3 content is observed in the northwestern part of the study area, where the aquifer is confined, the groundwater extraction rate is low and the cover layer is thick and composed of fine grain-size alluvium deposits.
The NO3 in the recharge zones is not monitored, because of the lack of monitoring wells in these zones.The nitrates' transport velocity in the unsaturated zone is determined by many factors, such as rock types, permeability, porosity and amount of groundwater recharge [129].The cover layer in the recharge zones of the study area is composed of coarse alluvium deposits of high porosity and permeability.Its thickness varies from zero to 5 m [49,53,54].As mentioned in Section 2.3, the hydraulic conductivity in the recharge zones of the aquifer varies from 180 to 300 m/day [57,58,60].Therefore, it is expected that the travel time of the nitrates in the recharge zones of the study area will be insignificant.
Figure 16 indicates that the variation trend of the nitrates content in the groundwater is in line with the variation trend of the vulnerability classes of the (final) vulnerability map.Once again, it can be said that the use of nitrates content as a validation technique supports the rational selection of SINTACS as a GWV model and the necessity of the performed statistical analysis, as well as a realistic approach concerning the effective weights assigned to the SINTACS parameters.

Groundwater Vulnerability Changes Over Time
The present map of the intrinsic GWV of the area under study may change over time.Possible changes may occur as a result of the addition of hydrogeological data and information, the addition of new monitoring wells, changes in land use, climate change, etc. Figure 16 shows that the higher the NO 3 content in the groundwater, the higher is the aquifer vulnerability and vice versa.The overall low level of NO 3 content in the whole aquifer could be due to the very low use of fertilizers in agriculture during the last 30 years.There are no studies on nitrates' travel time throughout the cover layer in the study area.The travel time of the nitrates depends on the water table depth, the hydraulic properties of the cover layer, the soil type, etc. [101].Bancheri et al. calculated that the average travel time of the nitrates through a clayey-rich cover layer and the overlaid soil in a case study in Italy was roughly 37 years [101].The higher the clayey content in the cover layer and soil, the longer is the travel time of the nitrates from the surface to the groundwater.Regarding the Tirana-Ishmi aquifer, the highest nitrates content is observed in the Tirana City area and its neighboring area, where the aquifer is phreatic, the water table is almost shallow, the cover layer is composed of medium grain-size alluvium deposits, and both the levels of the groundwater extraction rate and the amount of wastewater discharged into the local rivers' waters are high.The sources of nitrates in groundwater in this part of the study area are mostly sewage and industrial activities [29,100,103,128].The lowest level of NO 3 content is observed in the northwestern part of the study area, where the aquifer is confined, the groundwater extraction rate is low and the cover layer is thick and composed of fine grain-size alluvium deposits.
The NO 3 in the recharge zones is not monitored, because of the lack of monitoring wells in these zones.The nitrates' transport velocity in the unsaturated zone is determined by many factors, such as rock types, permeability, porosity and amount of groundwater recharge [129].The cover layer in the recharge zones of the study area is composed of coarse alluvium deposits of high porosity and permeability.Its thickness varies from zero to 5 m [49,53,54].As mentioned in Section 2.3, the hydraulic conductivity in the recharge zones of the aquifer varies from 180 to 300 m/day [57,58,60].Therefore, it is expected that the travel time of the nitrates in the recharge zones of the study area will be insignificant.
Figure 16 indicates that the variation trend of the nitrates content in the groundwater is in line with the variation trend of the vulnerability classes of the (final) vulnerability map.Once again, it can be said that the use of nitrates content as a validation technique supports the rational selection of SINTACS as a GWV model and the necessity of the performed statistical analysis, as well as a realistic approach concerning the effective weights assigned to the SINTACS parameters.

Groundwater Vulnerability Changes Over Time
The present map of the intrinsic GWV of the area under study may change over time.Possible changes may occur as a result of the addition of hydrogeological data and information, the addition of new monitoring wells, changes in land use, climate change, etc.
Currently, there are no monitoring wells within the recharge zones.Additional data from additional wells-in particular, monitoring wells within the recharge areas-can change the surface size and geometrical configuration of the areas of an extremely high vulnerability class.However, the vulnerability class of the recharge zones will remain extremely high.In general, due to additional data and information, the sub-zones with different vulnerability classes may change in surface size and horizontal geometrical configuration, but it is thought that the core of the prepared intrinsic vulnerability map will not change, because of the reliability of the used data and information and the careful selection of the assessment method.Given that the GIS-created model is a "living model", the use of additional data and information for each SINTACS parameter can improve the study outputs and reduce the hydrologists' subjectivism.
The study area consists mainly of urban and agricultural areas.There are no forests in the study area, the deforestation of which would imply changes in the soil media properties and infiltration [130].Population growth, intensive urbanization and increasing economic activity will intensify the extracted groundwater amount and therefore the risk to groundwater pollution.In the present GWV map, the sub-zone with high vulnerability corresponds to Tirana City and its surroundings, where the groundwater extraction rate is high and currently the municipal wastewater is discharged into the local rivers.The construction of a WWTP for Tirana City, as planned by the Tirana Municipality [27], will reduce this risk.But increases in the amount of groundwater extraction may affect the water table depth, which could slightly modify the present GWV map.However, it should be emphasized that in recent years, additional sources of drinking water for Tirana City are being sought outside the study area [43].
Climate change may affect the GWV [131,132] and therefore the outputs of the present study.Climate change projections for the Albanian Western Lowland [133], including the area under study, include precipitation and temperature rises.Changes in rainfall amounts and their seasonal distribution may impact the GWV.A temperature increase will increase evapotranspiration, which in turn may reduce the portion of surface water infiltrating to the groundwater.Therefore, it is expected that the projected climate change may affect the water table depth and infiltration.A statistical analysis helps to understand why and how climate change may affect the prepared vulnerability map.Table 9 indicates that small changes in water table depth and infiltration parameters have a great influence on the vulnerability index, due to the high variation coefficient.

Study's Limitations
This article has some limitations, which should be taken into consideration in future studies on the aquifer's vulnerability.The potential impact of densely urbanized areas on groundwater pollution could be considered as an additional parameter that should be taken into consideration through using a modified SINTACS method.Furthermore, as mentioned in the previous section, it is necessary to complete the groundwater monitoring network with additional monitoring wells, especially in the recharge areas, which correspond to the areas of an extremely high vulnerability class.

Conclusions
Based on the aquifer characteristics and the available data and information, SINTACS resulted in the most realistic groundwater assessment model for the Tirana-Ishmi alluvium aquifer.Validation through both a sensitivity analysis and nitrates content confirmed the rational selection of the SINTACS model and the reliability of the study's outputs and outcomes.The aquifer vulnerability to pollution assessment and the prepared vulnerability map indicate the following main conclusions: • The recharge areas show an "extremely high" vulnerability class (or degree) due to absence of both the unsaturated zone and soil media.This allows direct contact between the surface waters and the groundwater.

Figure 1 .
Figure 1.Hydrogeological map of the study area (according to the hydrogeological map of Albania, 1:200,000).

Figure 2 .
Figure 2. Geological map of the study area (according to geological map of Albania, 1:200,000).

Figure 2 .
Figure 2. Geological map of the study area (according to geological map of Albania, 1:200,000).

Figure 4 .
Figure 4.A 3D model of the aquifer.

Figure 4 .
Figure 4.A 3D model of the aquifer.

Figure 5 .
Figure 5. Map of the water table depth parameter.

Figure 5 .
Figure 5. Map of the water table depth parameter.

Figure 6 .
Figure 6.Map of the net recharge parameter.

Figure 6 .
Figure 6.Map of the net recharge parameter.

Figure 7 .
Figure 7. Map of the unsaturated zone parameter.

Figure 7 .
Figure 7. Map of the unsaturated zone parameter.

Figure 8 .
Figure 8. Map of the soil media parameter.

Figure 8 .
Figure 8. Map of the soil media parameter.

Figure 9 .
Figure 9. Map of the aquifer media parameter.

Figure 9 .
Figure 9. Map of the aquifer media parameter.

Figure 10 .
Figure 10.Map of the hydraulic conductivity parameter.

Figure 10 .
Figure 10.Map of the hydraulic conductivity parameter.

Figure 11 .
Figure 11.Map of topographic slope.(a) Topographic slope according to terrain map, (b) map of topographic slope parameter.

Figure 11 .
Figure 11.Map of topographic slope.(a) Topographic slope according to terrain map, (b) map of topographic slope parameter.

Figure 12 .
Figure 12.Preliminary GWV map of the Tirana-Ishmi aquifer at 1:50,000 scale based on theoretical weights of the SINTACS model parameters.

Figure 13 .
Figure 13.Mean effective and theoretical weights of the SINTACS parameters.

Figure 13 .
Figure 13.Mean effective and theoretical weights of the SINTACS parameters.

Figure 14 .
Figure 14.GWV map of the Tirana-Ishmi aquifer at 1:50,000 scale based on effective weights o SINTACS model parameters.

Figure 14 .
Figure 14.GWV map of the Tirana-Ishmi aquifer at 1:50,000 scale based on effective weights of the SINTACS model parameters.

Figure 15 .
Figure 15.GWV classes in the Tirana-Ishmi aquifer using theoretical weights of the SINTACS parameters, and the related nitrates content values.

Figure 15 .
Figure 15.GWV classes in the Tirana-Ishmi aquifer using theoretical weights of the SINTACS parameters, and the related nitrates content values.

Figure 16 32 Figure 16
Figure16shows the monitored nitrates content in the aquifer and the vulnerability classes after the model's submission to sensitivity analysis.Figures17-19show the vulnerability classes and the spatial and temporal variation in the nitrates content during the period 2015-2022[104].

Figure 16 .
Figure 16.GWV classes in the Tirana-Ishmi aquifer using effective weights of the SINTACS parameters, and the related nitrates content values.

Figure 17 .
Figure 17."Very high" GWV class and spatial and temporal distribution of the nitrates content.

Figure 18 .
Figure 18."High" GWV class and spatial and temporal distribution of nitrates content.

Figure 16 . 32 Figure 16
Figure 16.GWV classes in the Tirana-Ishmi aquifer using effective weights of the SINTACS parameters, and the related nitrates content values.

Figure 16 .
Figure 16.GWV classes in the Tirana-Ishmi aquifer using effective weights of the SINTACS parameters, and the related nitrates content values.

Figure 17 .
Figure 17."Very high" GWV class and spatial and temporal distribution of the nitrates content.

Figure 18 .
Figure 18."High" GWV class and spatial and temporal distribution of nitrates content.

Figure 17 . 32 Figure 16
Figure 17."Very high" GWV class and spatial and temporal distribution of the nitrates content.

Figure 16 .
Figure 16.GWV classes in the Tirana-Ishmi aquifer using effective weights of the SINTACS parameters, and the related nitrates content values.

Figure 17 .
Figure 17."Very high" GWV class and spatial and temporal distribution of the nitrates content.

Figure 18 .
Figure 18."High" GWV class and spatial and temporal distribution of nitrates content.Figure 18. "High" GWV class and spatial and temporal distribution of nitrates content.

Figure 18 .
Figure 18."High" GWV class and spatial and temporal distribution of nitrates content.Figure 18. "High" GWV class and spatial and temporal distribution of nitrates content.

Figure 19 .
Figure 19."Moderate" GWV class and spatial and temporal distribution of nitrates content.

Figure 19 .
Figure 19."Moderate" GWV class and spatial and temporal distribution of nitrates content.

Table 2 .
Groundwater vulnerability (GWV) classes in terms of the vulnerability index.

Table 3 .
Intervals of groundwater depth and SINTACS rating evaluation.

Table 4 .
Intervals of cover layer thickness, potential infiltration coefficient and SINTACS rate.

Table 4 .
Intervals of cover layer thickness, potential infiltration coefficient and SINTACS rate.

Table 5 .
Lithotypes of the unsaturated zone and the respective SINTACS range.

Table 5 .
Lithotypes of the unsaturated zone and the respective SINTACS range.

Table 6 .
Soil types in the study area and the respective SINTACS rates.

Table 6 .
Soil types in the study area and the respective SINTACS rates.

Table 7 .
Lithotypes of the saturated zone and the respective SINTACS rates.

Table 7 .
Lithotypes of the saturated zone and the respective SINTACS rates.

Table 8 .
Hydraulic conductivity and the respective SINTACS rates.

Table 9 .
Statistical analysis of SINTACS parameters.

Table 9 .
Statistical analysis of SINTACS parameters.

Table 10 .
GWV classes and the related surface areas for the Tirana-Ishmi aquifer according to the theoretical weights of the SINTACS parameters.

Table 10 .
GWV classes and the related surface areas for the Tirana-Ishmi aquifer according to theoretical weights of the SINTACS parameters.

Table 11 .
Map removal sensitivity analysis, removing one SINTACS parameter.

Table 12 and
Figure

Table 12 and
Figure

Table 13 .
GWV classes and the related aquifer's surface areas based on effective weights.

Table 13 .
GWV classes and the related aquifer's surface areas based on effective weights.