Determination of Groundwater Protection Zones of the Pozharan Wellfield Using Hydrogeological Modflow Model

This article concerns the calculation of the Groundwater Protection Zones (GWPZ) of the Pozharan/Požaranje wellfield. It shows the methodology of delineating water source protection zones with a hydrogeological computer model and serves as an example for further work in this field. The wellfield is located in the south-eastern part of Kosovo, about 1 km west of Viti/Vitia and is an important water supply source for the neighboring villages. The wellfield is located 35 km southeast of Prishtina/Priština, the capital city of Kosovo. A total of four public water production wells have been drilled into the aquifer for which the protection zones will be calculated. In order to delineate the Groundwater Protection Zones according to the Kosovar regulations, a groundwater model was set up to calculate the groundwater flow in the well field. Data has to be collected to create such a model. With help of previous studies and own investigations, the aquifer was identified. A large part of the work is finding observation wells (piezometers) in the study area and measuring its height and groundwater level. Afterwards, the model was calibrated. The model is capable of calculating flow paths and by means of particle tracking, it is possible to visualize where the water comes from. Adding the speed of groundwater flow, the time dependent zones can be drawn. Finally, the three protection zones were described as well as the proposed land use restrictions and the recommendations for land use planning were described. Several hazards to groundwater were identified and described inside those zones.


INTRODUCTION
Improved water source protection in piloted municipalities is one of the tasks of the RWSSP VI. Water source protection is a new issue in Kosovo and therefore a considerable challenge for water safety. The program has the task to define and delineate Groundwater Protection Zones (GWPZ) in pilot wellfield sites in Kosovo and to support all institutions and organizations to assume their role and tasks in the process to implement the protection measures.
In order to delineate groundwater protection zones 1, 2 and 3 according to the Kosovar regulations (MESP, 2017), a groundwater model was set up to calculate the groundwater flow in the well field. A so called 50-day-line was calculated to facilitate the implementation of a groundwater protection zone 2. For zone 3, the whole surface for the recharge of the well field was calculated.
The technical investigations require detailed surveys, starting from the geological, hydrogeological and morphological field work as well as monitoring of groundwater levels, realization of additional piezometers, calculation of hydrogeological parameters and several more. For the Pozharan well field, a hydrogeological computer model was established to understand the groundwater flow and to delineate the three necessary groundwater source protection zones. The technical understanding further requires a pollution survey and a water quality analysis.
The methodology focuses especially on achieving satisfactory result while modeling with Processing Modflow under limited data conditions. The authors aimed to describe how the water demand was estimated due to a lack of water production data and how the data from the existing private wells were used to understand the groundwater dynamics in the alluvial aquifer. The results show how the groundwater flow model leads to satisfactory results.
This article aimed to describe how these Groundwater Protection Zones are defined and designated based on a computer model as well as relevant geological and hydrogeological data.

General overview
The study area is located in the South-East of the Republic of Kosovo in the municipality of Viti/Vitia and is located between two rivers (see Figure 1).
The pilot project tackles the Pozharan/ Požaranje wellfield. The groundwater source in Pozharan/Požaranje stems from two wells which supply the village Pozharan/Požaranje with its 750 households and 3750 inhabitants and as well the village Çifllak/Čiflak with its 30 households and 150 inhabitants. This well field is also penetrated by two more wells which supply several more villages in the area (mainly Slatina/ Sllatinë). Figure 1 shows the geographic location of the study area including all relevant wells.
The study area is drained by Morava Binça, which collects all the small rivers. It has an average monthly flow rate of 6.7 m 3 /s. The surface of the basin is 1,560 m², whereas its average flow on the exit from Kosovo territory is 11.0 m³/s (PZHK, 2010). In Morava Binça, the maximum water flows occur in February, March and May, while the lowest are in August and September.
The surface water flows of this region, besides Morava Binça are the rivers Gërnçar, Debelldeh, Letnica, and Gelbush, etc. In the months with high water levels in February, March, May and during the months with maximum rainfalls, the rivers of Morava Binça and Glibusha flood the surrounding lands.

Meteorology
The explored region is characterized by an average continental climate, with relatively hot summers and moderately cold winters. The annual average rainfalls in Pozharan/Požaranje reach 820.6 mm, while the rainfalls in the surrounding mountains have to be considered higher. The months with highest rainfall are May, October

Extraction Rates of the wells
The extraction rate of the pumping wells is crucial when understanding the protection zone around the well. The extraction rate used for the model is based on the calculated optimum pumping rates of the wells. According to the operating water utility, the actual pumping rates are respecting these values. These numbers are as follows: The total water production capacity in the wells cumulates to 56 l/s (= 4838 m 3 /d). (Pruthi, 2008) (Pruthi, 2013) The Geology in the Study Area The wide region of the Pozharan/Požaranje aquifer consists of various lithologies of different ages -from the lower Paleozoic to Quaternary (see Figure 2).
The oldest rocks of the region are the different schists of the Lower Paleozoic. The crystalline schists build the basement of the ground part of the tectonic depression of Morava e Binçës (Anamorava) and the parts of its northern and southern sides. In the field area of the tectonic depression, the Paleozoic rocks are covered by the Tertiary and Quaternary sediments, while on the surface they appear on the hills of Karadak (Skopje) in the south and on the mountain hills of the Zhegovc in the north  Alluvial deposits (aQh) occur in all current water flows, especially in low flat regions and are represented by river sediments, sand and gravel deposits. They are widespread in the areas along the rivers Glibusha, Smira and Viti/Vitina.

The Hydrogeology in the Study Area
The wellfield of the Pozharan/Požaranje aquifers, is characterized as follows: Pliocene-Quaternary sediments (lake sediments, Proluvial (alluvial fan) and alluvial deposits) are predominant in the tectonic depression of Morava Binça. The central part of the plain of tectonic depression consists of alluvial deposits, river terraces. The quaternary sediments (Proluvial and alluvial deposits) represent an important intergranular aquifer.
The material of alluvial sediments, formed by the surface waters, is well rounded and crushed. These are mainly clastic sediments with high intergranular porosity, deposited as river terraces and in the river bed. They contain substantial groundwater reserves. The alluvial deposits are classified (by size) and are of heterogeneous composition. Alluvial sediments are unconfined aquifers with free groundwater surface. The alluvial aquifer has good hydraulic conductivity and effective porosity, recharged by infiltrations from the surface water. The groundwater level is relatively close to the ground surface. In some cases, leakages from other aquifers recharge of the considered aquifer, apart from the surface infiltration of rivers and rainwater. These waters are of poor mineralization and belong to the group of calcium-hydrocarbon waters.
The unconfined aquifer of alluvial sediments, accumulating considerable reserves of groundwater, is of great importance for the water supply; on the other hand, it is shallow and protected by a thin clay coverage, whereas waters of this aquifer are vulnerable to contamination. (KPMM, 2006)

Hydraulic Conductivity
On the basis of several pumping tests, the hydraulic conductivity of the area ranges between

Hydroisohypse
A total of 12 observation wells were identified and their water levels measured several times per year. These provide the groundwater levels and the directions of groundwater flow through the calculation of a hydroisohypse map. Figure 3 shows one example of the hydroisohypse map from summer 2018.

METHODOLOGY
An important part of the methodology to establish a water source protection zone is the field investigation and the understanding of the study area. These steps are standard hydrogeological procedures and shall not be explained in detail, since the focus of this article is the modeling. On the basis of the steps described above, the information presented in chapter 2 was discovered. It serves as a basis in setting up the computer model. The chosen software was MODFLOW developed by the US Geological Survey (Harbaugh, 2005). MODFLOW is commonly used model in delineating Wellfield Protection Zones (Forster et al 1997). It is considered an international standard for simulating and predicting the groundwater conditions and groundwater/surface-water interactions. The Graphical User Interphase (GUI) is Processing Modflow (Chiang & Kinzelbach, 1998). In addition, PMPATH was used, which is running transport models independently of MODFLOW. It retrieves the groundwater simulations results from the MODFLOW software to calculate the groundwater paths and travel times (Chiang & Kinzelbach, 2003). The MODFLOW model is an efficient tool for determination of the water protection zone, which is essential in water resource management (Goodarzi et al., 2019). It is user-friendly software popular among the hydrogeologists for its features (Hariharan V, 2017).

Input Parameters
The quality of the input data in groundwater models has a significant effect on the model results

Grid
The mesh extends 4,450 m from west to east and 4,200 m from south to north. The grid is more detailed around the wells and more rough in less important areas of the model. The cells in the wellfield area are scaled down to 10 m, while the cells far away have a size of 50 m. In total, the model has 127 rows and 118 columns, which amounts to 14,986 cells.

Layer
In order to define the model thickness, each layer is determined by Top of Layers (TOP) and Bottom of Layers (BOT) files. For this purpose, model 3 layers are used. Layer 1 embodies the topsoil and Layer 2 and Layer 3 comprises the aquifer. For each layer, an elevation grid for TOP and BOT is created. BOT of Layer 1 (BOT 1) is the same as TOP of Layer 2 (TOP 2), so 4 elevation grids are needed. TOP of Layer 1 (TOP 1) represents the surface level. TOP 2 (or BOT 1) is the lower level of the soil. Each grid is based on different data sets. The data gathered during the measurement campaign is used for the ground level. During prior investigations of the well field, the wells were topographically surveyed and integrated to the grid. Each point has X, Y and Z values which are entered to the Surfer software to create a grid through the Kriging interpolation method. The grid is compared with a topographic map to find implausible levels of the grid and correct them.

Time
Limited data availability for the meteorological and hydrological events requires the steady state simulation. This steady state simulation uses average annual data for precipitation and water extraction from the pumps. It is not possible to calculate the extreme events like a flooding event.
The length of stress periods and time steps is not relevant to steady state simulations, only for transport. However, it is important to stick to the same unit, especially for parameters like hydraulic conductivity.

Initial hydraulic heads
Modflow software requires the definition of hydraulic heads as starting values in steady state simulations. The initial heads are plotted in the constant head cells and represent the actually measured values from the field surveys.

Horizontal hydraulic conductivity and transmissivity
The horizontal hydraulic conductivity is a value describing the permeability of the underground. Hydraulic conductivity calculation is based on pumping tests. However, the geologist has to compare the values with standard values for the various types of rock underground and check for plausibility.

Effective Porosity
Porosity describes the ratio of open volume and solid volume of rocks and soils. The effective porosity is less than total porosity, as parts of the fluid in pore space are immobile. There are no laboratory values available so far, so that the effective porosity has to be estimated. The literature values provide a wide range of values for each type of rock. The present model is calculated with 0.1 of effective porosity, as described, because of the arithmetic mean for gravelly and sandy soil.

Boundary Condition
In Modflow, each cell requires an IBOUND array that defines whether a cell is active, inactive or provides a constant head. Active cells are used in the aquifer area of the model. The groundwater levels and flow are calculated in such cells. In turn, inactive cells are used outside the aquifer area, where no groundwater flow is expected.
The initial hydraulic head values remain fixed throughout the simulation in constant head cells. In this model, a constant head boundary exists in the areas outside the influence of the well field. The model boundaries are either inactive or constant head cells. In order to simulate the inflow and outflow of the model area active parts of model boundaries are set to a constant head.

Packages
Packages are modules to the basic Modflow software providing the options to simulate aspects of the whole water cycle.

Recharge Package
The Recharge package simulates the distributed recharge to the groundwater system. The recharge is assigned to each vertical column of cells. The water added to the system through precipitation is simulated with the recharge package. It is estimated that about 36% of the precipitation reaches the aquifer and acts as recharge. The annual precipitation in 2017 amounted to 820 mm; therefore, an average daily recharge of 8·10 -5 m was estimated.

River Package
The Morava river and Gilbusha river run across the study area and need to be included into the model. The package simulates the flow between the river and the groundwater system. The input parameters are hydraulic leakage coefficient of river bed, head of river, elevation of the riverbed bottom, width of the river.

Well Package
The well package is used to simulate the water abstractions by wells. As detailed data is not available due to the absence of water meters on well heads, an average over the annual production of the two abstraction wells was used. The total annual abstraction was given as approx. 215.350 m³, resulting in a daily abstraction of 295 m³/d by each of both wells.

Calibration
Calibration means adjusting the model parameters in an allowed range until the modeled results correspond realistically to the measured results. The soil parameters and weather dependent data are given and should not be modified significantly. Figure 1 compares the final groundwater levels of the model to the measured levels. The level is in average 0.4 m above the measured level. The reason most probably is a decreased precipitation rate in the summer of 2018.
The production of the model was not changed to reduce the risks for negative misinterpretation. A dynamic behavior of the observation wells due to specific condition like temporary abstraction, irrigation activities etc. further contributes to the discrepancy between modeled and measured average values.
Another point of uncertainty is the unknown daily variation of production, providing direct changes in the surrounding wells.

RESULTS AND DISCUSSION
The main result from modeling is the different water protection zones. Zone 1 is the direct protection of the well within 10 m around the well and does not need to be modeled. Zone 2 refers to the line from which the water needs 50 days travel time to reach to the well. It is thus called the 50-day line, which is the first result. Figure 5 provides the resulting GWPZ 2 enveloping curve printed on the topographic map of the study area.
Two circular lines around each of the production wells are formed. Each of the two circles has an approximate 70 m diameter which represents the minimum extension of the GWPZ 2.
The Recharge area of the well field corresponds to GWPZ 3. The recharge area is the surface area from where flow paths of infiltrating surface water (mainly precipitation) turn towards the abstraction wells and reach there in any undefined time period ( Figure 6). In reverse conclusion, the infiltrating water of the surface area outside this recharge area will not travel towards the abstraction well, but will be transported away from the abstraction wells.
The Recharge has an 8-shape and extends approx. 2.5 km in the North-South direction. The southern East-West extension is 2 km long, the northern E-W extension is around 1.5 km. The northern end of the recharge area already includes the most southern houses of Pozharan/ Požaranje village. The remaining area of GWPZ 3 is an agricultural area.

CONCLUSIONS AND OUTLOOK
Groundwater models are useful tools of water resource planning and management (Uddameri et al. 2017). As a general result of this modeling of the well field, it has to be concluded that two GWPZ 2 have to be declared. Each of them with approximately 70 m of diameter around the production wells. The GWPZ 3 is a larger area including the southern houses of the Pozharan/ Požaranje village. In this case, one area is the recharge area for all four wells. This hydrogeological investigation is the first step in implementing the Groundwater Source Protection zones. The results from this analysis have been included in a technical report (RWSSP-VI, 2019) in which the 50-day line and the recharge area are presented.
It is the responsibility of the Ministry of Environment and Spatial Planning of Kosovo to implement the Groundwater Source Protection Zones based on the calculated areas. For each zone, different measures are to be taken in the future. These measures are: • GWPZ 1: The wells need to be fenced according to regulations. • GWPZ 2: The protected area has to be declared limited and agricultural activities shall be restricted according to the regulations in force. Contaminated water, waste and other hazards to the groundwater can infiltrate into the aquifer easily through holes (rural wells or piezometers). It is necessary to close these or protect these against infiltration, especially if they are close to the wells. • GWPZ 3: All over the river banks of Glibusha river and along streets and pathways, residents dump their waste. Waste of different varieties has been observed: waste of house construction and other infrastructure, domestic waste, cans of motor oil and paint etc. Such contamination cannot be tolerated in the area of GWPZ 3. Besides the slow leakage of pollutants into the groundwater, the waste floats around during the flooding periods. Therefore, these floods have to be considered as potentially dangerously contaminated surface water. The direct Figure 6. Recharge area - GWPZ 3 infiltration via the open rural wells and piezometer drillings has to be considered a high risk to water quality. All waste dumps inside the GWPZ 3 have to be removed and cleaned.