Delineation of Well Head Protection Areas for the Public Wells in the Ferizaj Region (Kosovo) with Limited Data Availability ()
1. Introduction
In the past decade, the Rural Water and Sanitation Support Programme in Kosovo (RWSSP) funded by the Swiss Agency for Development and Cooperation (SDC) and implemented by Dorsch International Consultants GmbH and Community Development Initiatives (CDI) has considerably increased the access to public water supply of the rural population of Kosovo. However, groundwater protection is still a challenge [1] [2]. While protecting wells with modeled groundwater protection zones is similar to German standards [3] [4] [5] and is well described in Kosovar regulations [6]. The full implementation of modeling-based groundwater protection areas, however, is yet to be achieved [7]. The programme is thus supporting the Regional Water Companies of Kosovo to protect groundwater resources used for public water supply.
This study describes the development of a consistent hydrogeological model and the delineation of well head protection areas for the public wells in the Ferizaj region by applying a mathematical groundwater flow model under limited data availability. However, data availability with a high degree of accuracy is not always available which, at first sight may limit the use of mathematical groundwater models. The study shows that even when not all data are available, it is possible and necessary to use mathematical groundwater models to delineate well head protection zones.
Data collection, field work, and evaluation of the data were carried out by Dorsch International Consultants GmbH and Community Development Initiatives (CDI) in close cooperation with the Regional Water Company Bifurkacioni in Ferizaj in 2018/2019.
2. Material and Methods
2.1. Study Area
The study area is in the south of Kosovo in the region of Ferizaj, which includes the municipality of Ferizaj (Figure 1). Within the region, there are numerous small villages in the valley. The study area comprises the upper part of the catchment area of the Nerodime river including the alluvial valley plain and the surrounding mountains. The elevation of the valley plain ranges from south to north from 530 to 600 m a.s.l., the average width of the valley is 5 km. The valley is surrounded by mountains with a maximum elevation up to 1500 m a.s.l. The study area has a total area of 195 km2 with a maximal width of 20 km and maximal length of 10 km.
In the study area, 14 drilled wells have been implemented. Six wells are used as production wells for drinking water supply. Of these, two wells each are located near the villages Varosh, Gerlice and Begrace.
2.2. Climatic and Hydrologic Conditions
The region is characterized by an average continental climate, with relatively hot summers and moderately cold winters. Based on meteorological/climatic data, the region of Ferizaj is classified as climate “C” according to the nomenclature defined by Köppen and Geiger [8].
![]()
Figure 1. (a): Location of the study area. (b): Topographic map of the study area. (c): Map of the study area showing the location of Ferizaj city, villages, Varosh-, Gerlice-, Begrace-well field, observation wells and the extent of the alluvial aquifer and the adjacent hard rocks.
The mean long-term annual precipitation for the time period 2002 to 2018 is P = 737 mm/a, with 357 mm for the hydrological winter months (October to March), and 380 mm for the summer months (April to September). February is the driest month with a precipitation of 36 mm/month, May the wettest month with a precipitation of 86 mm/month [9]. The mean annual temperature is T = 11.0˚C, the mean winter temperature is 4˚C, and the mean summer temperature is 17.3˚C [9].
The average long-term annual runoff measured in the Nerodime river at a gauge station some kilometer south of the study area is Q = 1.2 m3/s (time period 2015 to 2018). The corresponding catchment area is 214 km2, so the total runoff is calculated at R = 177 mm/a (=5.6 l/s∙km2) [9]. The mean annual evapotranspiration is thus calculated from the water budget at E = 560 mm/a. The mean long-term water balance for the study area is:
P = E + R
737 = 560 + 177 [mm/a]
23.4 = 17.8 + 5.6 [l/s∙km2]
100 = 76 + 24 [%]
The total runoff R is made up of surface runoff and groundwater runoff. The groundwater runoff corresponds to the groundwater recharge Ig. The groundwater recharge was calculated from the discharge data using the “MoMNQ-Verfahren” by Wundt [10]. In this method, the average groundwater outflow is determined by using the average of minimum monthly streamflow. For the study area groundwater recharge was calculated to Ig = 122 mm/a (=4 l/s∙km2).
2.3. Geological and Hydrogeological Settings
The study area is part of the so-called “External Vardar Subzone” bordering the “Central Vardar Zone” to the east [11] [12] [13]. The main plains containing the studied aquifer are part of the “Kosovo Basin”.
The emergence of this intra-mountainous basin started in Neogene. The deposits consist mainly of unconsolidated fluvial sediments of Pleistocene and Holocene age formed in a floodplain environment of a meandering river. The surrounding hard rocks comprise of consolidated sediments (sand-, lime-, clay-, marl-, and siltstone), magmatic, and metamorphic rocks from Palaeozoic to Tertiary age. The region is tectonically stressed. The main tectonic faults strike parallel to the axis of the valley. The hydrogeological map of the study area is shown in Figure 2 [14]. Figure 3 shows a geologic cross-section of the alluvial aquifer.
The alluvial aquifer can be characterized as an unconsolidated porous aquifer [15]. The aquifer consists of an alteration of gravel, sand, silt, and clay. Its lateral extent from west to east ranges from about 3 to 7 km. The aquifer thickness varies between 160 to 260 m. The average thickness is about 200 m. The depth to groundwater is only a few meters. The aquifer is underlain and surrounded by hard rocks which form the bottom of the aquifer. The aquifer is in the upper part in an unconfined state while deeper parts of the aquifer are considered to be in a confined state.
The adjacent and underlying hard rocks can be characterized as aquitards. The hydraulic conductivity and storativity of these hard rocks is very low compared to the alluvial sediments even though the karstic parts in the limestones may be considered of good permeability.
Hydraulic conductivity: The estimation of the hydraulic conductivity K was done by the analysis and evaluation of 14 pumping tests carried out in the study area. The names of the wells and the estimated K-values are shown in Table 1. The K-values range from K = 7.52 × 10−6 to 3.30 × 10−5 m/s. The average K-value computed as the geometric mean is K = 2.00 × 10−5 m/s.
Figure 4 shows a map of the spatial distribution of the hydraulic conductivity values in the alluvial aquifer and the corresponding empirical frequency distribution. The histogram shows that the logarithms of the K-values are normally distributed [16], which was approved by the Kolmogorov-Smirnov Test (D = 0.22451, p = 0.05407). The standard deviation and the range are comparatively low, which indicates a homogeneous structure of the aquifer.
![]()
Figure 2. Hydrogeological map of the study area with location of the cross-section through the alluvial valley sediments. Modified excerpt from the hydrogeological map of Kosovo [14].
![]()
Figure 3. Hydrogeologic cross-section of the alluvial aquifer (line of cross-sections is shown in Figure 2. Note: the vertical scale is exaggerated).
![]()
Figure 4. Map showing the spatial distribution of the hydraulic conductivity (K) in the alluvial aquifer and the corresponding empirical frequency distribution (
= mean, σ = standard deviation, n = number of samples).
![]()
Table 1. Hydraulic conductivity (K) from pumping test analysis for the wells in the alluvial aquifer.
Transmissivity: The corresponding transmissivities estimated from pumping test data are in the range from T = 3.5 × 10−4 to 2.4 × 10−3 m2/s. The average T-value is T = 9.8 × 10−4 m2/s.
Porosity: Within the framework of this study no special investigations for the estimation of total and effective porosity were carried out. However, it can be assumed that the average total porosity of the alluvial sediments is in the range of n = 10% to 20% and that there is no significant difference to the effective porosity nf.
Groundwater dynamics: The main groundwater flow in the alluvial aquifer is directed towards the Nerodime river, which represents the receiving river during average flow conditions. A significant inflow from the river into the alluvial aquifer occurs only during flood events. The values of the groundwater heads ranges from about 600 m a.s.l. in the north to 520 m a.s.l. in the south. The hydraulic gradient under natural flow conditions is in the range of I = 5 × 10−3 to 1 × 10−2. The natural groundwater flow velocity ranges from u = 0.1 to 0.5 m/d. Besides groundwater recharge through precipitation and infiltration in the outcrop area of the alluvial aquifer significant additional recharge occurs by groundwater inflow from the slopes of the adjacent mountains. The total groundwater volume in storage can be estimated from the aquifer volume which is about 11 km3 by multiplying this volume with the porosity of the aquifer. Assuming a porosity of n = 0.1 to 0.2 the total groundwater volume in storage within the alluvial aquifer is in the range of 1.1 to 2.2 km3.
Groundwater abstraction: Within the study area actually three well fields for the public water supply exist. These well fields are from north to south the Varosh-, Gerlice- and Begrace-well field. Each well field consists of two pumping wells. Information about the wells is given in Table 2.
The wells extract groundwater from the alluvial aquifer in the valley. The estimated actual groundwater abstractions are Q = 7.5 l/s for each well. The total groundwater abstraction for public water supply is Q = 45 l/s. Besides the public wells numerous private wells exist. These are shallow dug wells which are mainly
![]()
Table 2. Name, location, ground level, groundwater level, well depth and pumping rate for the public wells in the study area.
used for private irrigation of the gardens around the houses. The total groundwater extraction through these wells is estimated to be low and therefore neglected in the water budget analysis.
Observation wells: Since 2018, the RWSSP started the investigation of the Ferizaj well fields. During this campaign a total of 80 private wells and piezometers were identified, mapped, measured, and monitored. The main goal was to determine the groundwater level at these observation wells to achieve a reliable basis for the construction of the groundwater head contours. The location of the observation wells is shown in the maps.
3. Mathematical Groundwater Flow
3.1. Choice of Mathematical Model and Its Boundary Conditions
For the delineation of the well head protection areas and the quantification of the water budget a 2-dimensional steady-state groundwater flow model assuming average flow conditions was designed. As the well drawdown in the alluvial aquifer is small compared to the total thickness of the aquifer confined aquifer conditions were assumed.
For the numerical simulation of groundwater flow the program Modflow [17] [18] was used. Pre- and post-processing was done with Processing Modflow Version 8.047 [19] [20]. The delineation of the catchment areas of the wells and the corresponding 50-day isochrones of the production wells were calculated with the advective transport model PMPATH [21].
Model area: The model area comprises the upper part of the catchment area of the Nerodime river including the alluvial aquifer and the hard rocks of the adjacent mountains (Figure 5). That means that the model area includes the total catchment area of the production wells. This approach guarantees a feasible estimation of the total groundwater budget, especially the estimation and distribution of the groundwater inflow from the adjacent hard rocks of the mountain area into the alluvial aquifer. The actual model area is limited to the extent of the alluvial aquifer.
Model grid: A finite difference grid with 131,898 square cells with a constant cell size of 50 × 50 m over an area of 329.745 km2 was selected for the flow model. The area of the active cells is 190.273 km2. There are 494 cells in the x-direction (west to east) and 267 in the y-direction (north to south). This discretization gives a clear-cut idea of the groundwater flow field in the entire model area.
Boundary conditions: The boundary conditions are orientated on the natural hydraulic boundaries (Figure 5). The watershed of the Nerodime river is simulated as a no flow boundary. The southern boundary within the alluvial is a fixed head boundary set to 523 m a.s.l. following the observed groundwater head contour line for average groundwater flow conditions at this location. The river reaches within the alluvial are simulated by a leakage boundary, while the rivers, creeks and springs in the hard rock area are simulated as drains. This means that
![]()
Figure 5. Map of the model area showing the boundary conditions for the flow model, the location of the well fields, the location of the observation wells, the horizontal extent of the alluvial aquifer, and of the hard rocks.
in the first case effluent and influent conditions can be simulated, while in the second case only an outflow of groundwater into the rivers respectively into the drains is possible.
3.2. Calibration of the Groundwater Flow Model
Due to the limited data availability a simple calibration strategy was applied. Within the model area two homogenous zones were defined. One zone corresponds to the alluvial valley, the other zone to the area comprised by hard rocks. For each zone the transmissivity of the aquifer and the groundwater recharge from precipitation was estimated. While the data for the alluvial aquifer are from field investigations, the values for the hard rocks are plausible estimates. Table 3 shows the average and the range of the model parameters transmissivity and groundwater recharge from precipitation for the two zones.
The total abstraction rate of the 6 public wells is Q = 4.5 × 10−2 m3/s (=7.5 × 10−3 m3/s per well). The leakage factor respectively the hydraulic conductance of the river bed and drains were kept constant during the calibration process. For the hydraulic conductance of the river bed a value of Criv = 5 × 10−4 m2/s was estimated, and for the hydraulic conductance of the drains a value of Cdrain = 2 × 10−1 m2/s. The high value of drain hydraulic conductance guarantees that the drain acts as a perfect drain.
![]()
Table 3. Input parameters for the calibration process (average, min., max. value): Transmissivity (T), groundwater recharge from precipitation (Ig) for the two zones (alluvial aquifer, hard rocks).
During the calibration process transmissivity and groundwater recharge were varied within the defined ranges by trial and error (Table 3). The simulation results were compared with measured head data from observation wells in the alluvial sediments. Observed head data from the hard rocks were not available, so that a classical calibration of the groundwater heads in this zone was not possible. Constraints for modelling were that the simulated heads are smaller than ground level, the hydraulic gradient is much smaller than the topographic gradient, and the inflows into the alluvial aquifer are plausible. However, the simulated head distribution in the hard rock area is not unrealistic. It is controlled by the drainage network.
The values for the best parameter fit from the calibration and sensitivity analysis are displayed in Table 3 (column: best fit calibration). The comparison shows that there is no big difference between the measured field data and the calibrated data for the transmissivity and groundwater recharge in the alluvial aquifer. The calibrated values correspond nearly the average values from the field data.
Figure 6 shows the corresponding scatter plot (hcalculated versus hobserved). It can be seen that the deviation between the calculated and observed heads is reasonable taking into account the simple zonation approach of the model in only two hydrogeological units in respect to transmissivity and groundwater recharge.
The calculated mean square deviation of the fit is MSD = 27.62 m2 (the MSD is defined as the mean squared error between observed and calculated heads). The error can be considered compatible with the aims to be pursued in this research.
Figure 7 shows the calibrated piezometric head distribution in the alluvial aquifer for steady state conditions. It can be seen that the groundwater dynamic within the aquifer is strongly dependent from the inflow from the hard rocks and the outflow in the Nerodime river which acts as gaining river, that means that the biggest part of groundwater discharges into the river.
The water budget for steady state flow conditions is shown in Table 4 and in Figure 8. The arrows correspond to the inflows and outflows. Inflows are counted positive, outflows negative.
![]()
Figure 6. Scatter plot: comparison of calculated and observed piezometric heads.
![]()
Figure 7. Simulated piezometric head distribution and general groundwater flow direction in the alluvial aquifer for steady state flow conditions.
![]()
Figure 8. Schematic sketch of the water budget for steady state flow conditions for the alluvial aquifer.
![]()
Table 4. Groundwater budget for the alluvial aquifer for steady state conditions.
From the total volume of groundwater in storage V and the total flux in or out of the alluvial aquifer Q the mean residence time (Tr) of the groundwater in the aquifer can be calculated to Tr = 132 to 263 years.
4. Results and Discussion
On the basis of the calibrated groundwater flow field the groundwater protection areas of the wells were calculated. For the calculation of the 50 days isochrone and the travel times an effective porosity of nf = 10% and an average aquifer thickness of m = 200 m was supposed.
For the delineation of the groundwater protection zones, the licensed average groundwater abstraction rates for the pumping wells are required [3]. In the present case, however no reliable data on groundwater extraction were available. Therefore, an estimation of the expected pumping rates was achieved by applying standard procedures for a drinking water network. By calculating the water demand until 2034 for the supply area and taking into consideration possible water losses the total abstraction rate of the 3 well fields was estimated to Q = 45 l/s. Finally, after checking the specific capacity and pump size of the wells the discharge rate of the six pumping wells was set to Q = 7.5 l/s for each well.
The simulated extent of the well head protection areas for the well fields within the alluvial aquifer on the basis of the simulated steady state flow field is shown in Figure 9. In accordance with the general groundwater flow direction, the catchment areas run from west to east. The width of the catchment areas ranges from about 1.2 to 3.4 km. The groundwater travel time from the boundary of the alluvial aquifer to the well fields is in the magnitude of 25 to 40 years.
It can be seen that all well fields get a significant boundary inflow from the adjacent hard rocks. This inflow rate is of the order of about a third to half of the pumping rate of the wells. This means that part of the groundwater extracted by
![]()
Figure 9. Extent of the groundwater catchment areas for the well fields within the alluvial aquifer simulated on the basis of the steady state flow field. The red arrows show the boundary inflow rates from the adjacent hard rocks.
![]()
Figure 10. Example of groundwater protection area with 50-day isochrones for Gerlice well field.
the wells is replenished in the outcrop area of the hard rocks and not within the catchment area of the alluvial aquifer. Due to missing data, the delineation of the catchment areas within the hard rock area is not possible at the moment. A simple delineation only based on the groundwater recharge rate and the corresponding recharge area in the outcrop of the hard rocks is too inaccurate. Therefore, additional investigations of the hydrogeological conditions in the hard rocks are recommended.
Besides the simulation of the catchment areas, the extent of corresponding zone II of the well head protection area was calculated. The groundwater protection area zone II is for the protection against bacteriological pollution. The size of this area is defined by a 50-day isochrone around the pumping well. The groundwater travel time from the border of this zone to the well is 50 days. As an example, the 50-day isochrones and the corresponding pathlines for the two wells of the Gerlice well field is shown in Figure 10. The area of each zone is approximately 1600 m2, the maximum length and width are 44 and 50 m.
5. Conclusions
The study showed that for the technical delineation of well head protection areas the application of groundwater models is strongly recommended. A consistent quantitative interpretation of all hydrogeological data even under limited data availability is only possible with a mathematical groundwater model.
Furthermore, the analysis of the data with a model shows shortcomings in the data basis as well as in the understanding of the hydrogeological model and helps in the planning of additional investigations.
For the state of Kosovo where the protection of the groundwater resources is currently being carried out groundwater protection areas based on groundwater models will lead to scientifically proven protection areas which will have a higher public acceptance.
Acknowledgements
This study was kindly supported by the “Rural Water and Sanitation Support Program VI in Kosovo (RWSSP)” funded by the Swiss Agency for Development and Cooperation (SDC) and implemented by Dorsch International Consultants GmbH and Community Development Initiatives. The authors thank SDC for the possibility to develop this project and write this article. The article does not reflect the views of the above institutions.
Further, we want to thank the project team for their support and guidance during this study above all Thomas de Beyer, Hajrije Morina and Alan Brown. We further express our gratitude to Valdet Syla and Jeton Shabani from RWC Bifurkacioni for their support, contribution with data and valuable discussions.
During the programme a full technical report on the water source protection zone was developed for the Ferizaj area. For academic interest, some calculations of this report were revisited. This publication therefore serves as interesting additional information. For legal implementation of the zones, however, only the information in the technical report and not the ones from this article are relevant.