Determining the Groundwater Balance and Radius of Influence Using Hydrodynamic Modeling : Case Study of the Groundwater Source Šumice in Serbia

A groundwater flow model was developed to simulate groundwater extraction from the public water supply source of the City of Kikinda. The hydrodynamic model includes the municipal groundwater source of Kikinda (Šumice and the Jezero Well), but also an extended area where there are groundwater sources that provide water supply to three factories: (MSK, TM and LŽT Kikinda). Hydrodynamic modeling, based on the numerical method of finite differences will show the groundwater balance of the sources in the extended area of Kikinda. The impact of the industrial water sources on the regime of the public water supply source will also be assessed. The radius of influence of the groundwater source is determined by simulating the travel of conservative particles over a period of 200 days.


INTRODUCTION
Kikinda Municipality is located in the northeastern part of the Province of Vojvodina and occupies a land area of 782 km 2 .According to the most recent census (2002), the population count is 41,935 and the population growth rate -5.7%.Over the past several years, geological and hydrogeological research has been conducted in the extended area of Kikinda's groundwater source, aimed at providing the required amount of quality groundwater to meet the water demand of the population.The study area is situated in the southern part of the Pannonian Basin, which is a lowland that occupies the north of Serbia (Figure 1).The City of Kikinda extends over an area of 13.5 km 2 in the spacious Banat Plain, where the highest elevation is 83 m and the lowest 76 m.Water supply is provided to Kikinda by means of 11 wells.The tapped water-bearing horizon is at a depth of about 250 m, formed in Quaternary sands whose thickness is about 50 m.In the vicinity there are also groundwater sources for industrial water supply, which have a certain influence on the regime of the public water supply source.
Nowadays, water supply is a highly complex and significant issue in every society.It is of vital importance to control and monitor the operation of each water source.Unfortunately, there are many problems.Even when the required legislation is in place, it is sometimes impossible to implement it in practice, not only in the case of existing sources that have been in service for years, but also when it comes to opening new sources.The best research approach includes hydrodynamic analysis in the early stages of planning and sound monitoring in all stages.Many authors have recently focused on this subject from different perspectives: the hydrodynamic aspect of hydrogeological research [1][2][3], water source optimization and sustainable development of water supply systems [4][5][6], drinking water supply [7], different types of factors and processes that influence water source performance [8], water quality [9,10], and protection of water sources [6,11].
The groundwater balance of the groundwater sources was determined by hydrodynamic modeling.The impact of industrial groundwater sources on the municipal groundwater source was also assessed.The particle tracking method [12] was used to determine the radius of influence of the groundwater sources over a period of 200 days.The 200-day period reflected the travel time to the wells and represented the third and widest sanitary protection zone [13].Particle tracking can be used for various purposes.The backward particle tracking method with an uncertainty analysis of porosity, applying a Monte Carlo approach, with Geographic Information System (GIS) support, is a useful tool for delineating groundwater protection zones [14].Particle tracking simulations can be used to determine travel times from recharge to discharge areas along identified flowpaths [15].Using a standard numerical flow and transport code and a technique based on adjoint theory, and combining forward-in time and backward-in time transport modeling, it is possible to determine the impact of potential contaminant sources at unknown locations within a well capture zone, including the expected times of arrival of a contaminant, the dispersion-related reduction in concentration, the time taken to breach a certain quality objective and the corresponding exposure times [16].

GEOLOGICAL SETTING OF THE STUDY AREA
Extensive geological research (exploratory drilling, geophysical activities, paleontological investigations) has resulted in a relatively robust geological database [17].The lithological units comprise Pliocene and Quaternary sediments.The base of the Quaternary sediments is made up of shallow lacustrine sediments of the Upper Pliocene.They are built of sandy siltstone, gravelly and sandy siltstone, dark gray siltstone and organogenic silty sands.Carbonaceous interbeds and organic matter are regularly present in these sediments.The color of the sediments ranges from dark gray to light gray and olive gray.Pleistocene and Holocene deposits were abstracted during the Quaternary period.The earliest anthropogenic formations were built from Lower Pleistocene riverine and lacustrine sediments.Riverine/bog sediments of the Mindel period constitute an immediate sequel, whose continuity extended during the Mindel-Riss interglacial period.The Upper Pleistocene is comprised of riverine/bog sediments of the riverbed facies and floodplain facies, which constitute the Varoška Terrace deposits of the Würm age.Sedimentological analyses have shown that the riverine and lacustrine sediments were made of sands, silty sands and alevrites (sandy, rarely clayey and gravelly), while the riverine/bog sediments were composed of clay alevrites and alevrite sands with gravel.
Additional units were abstracted in the post-glacial period, during the Holocene: abandoned channel facies, floodplain facies, bog sediments, alluvia, and the youngest unit deposited since the Holocenethe beach facies of the Tisa River.Sedimentological analyses have shown that they were mostly made up of alevrites, sand, gravel and loess.

HYDROGEOLOGICAL OVERVIEW
A large number of wells (Figure 2) have been drilled in the study area [18].According to the type of porosity and hydrodynamic characteristics, there are two types of aquifers in this region: confined and unconfined, both in Quaternary sands.
The unconfined aquifer was formed in Quaternary alluvial sands.The depth of the sediments is about 30 m.The aquifer is rechared by infiltration of precipitation and surface water, since there is a good hydraulic contact between the aquifer and the river.The aquifer is drained naturally, at times of low river stages when groundwater flow is directed towards rivers, and artificially -through drilled and dug wells.
The confined aquifer was formed within several water-bearing complexes of Quaternary sands.Although this aquifer comprises a number of water-bearing layers, it was not possible to separate them on the basis of geometry, hydraulic interactions and physico chemical characteristics of the water.The reason for this is the evolution of the terrain, or, more precisely, there are lacustrine, bog and riverine sediments for which there are no reliable data about horizontal and vertical distributions and interbeds of semi-permeable and impermeable strata.An important characteristic of these layers is the constant alternation of sands, slightly gravelly sands, with loessial clays, sandy clays, and clays (Figure 3). Figure 3 shows a characteristic lateral alternation and thinning of impermeable layers and the water-bearing layer.The layers of sand, which are in places more than 30 m thick (although typically less than 10 m), represent the basis of certain polygenic packages.Overlying the sands, there are microgranular sediments: loamy sand, loessial clays, clays, and sandy clays.The transitions to fine-grain sediments are often gradual.Based on the grain-size distribution of the water-bearing horizons, which also include the studied aquifer, the material of the water-bearing formation is represented by fineand medium-grain sands (uniformity coefficient: Cu = d60/d10 = 2-3 mm with d50 = 0.1-0.3mm, where d is particle/grain size: d10, d50 and d60 represents a grain diameter for which 10%, 50% and 60% of the sample will be finer than it or using another words, 10%, 50% and 60% of the sample by weight is smaller than diameter d10, d50 and d60).According to pumping test data, the transmissibility coefficient in the extended area of Kikinda, depending on the locality in question, is in the range of (2.5-10) × 10 -2 m 2 /s.
Table 1 shows the hydrogeological parameters previously determined for the Kikinda area [19].The aquifer extends beyond the study area.The aquifer is probably recharged on the slopes of the Carpathian Mountains, the Vršac Mountains and Mt.Fruška Gora (east of the study area).It is also recharged from surface water in the southwestern part of the area, where the depth of the layers which from the aquifer is relatively small and allows a hydraulic contact between the river and the aquifer over the alluvial layers.The main types of aquifer drainage are: artificial drainage (through groundwater extraction) and leakage into overlying semi-permeable strata (given that the aquifer is confined), and from thereto the top aquifer.In lithological terms, aquitards are generally represented by clays and sandy clays (Figure 3).

CONCEPTUAL (HYDROGEOLOGICAL) MODEL
The conceptual model of the hydrogeological system is based on geological information on the boreholes (Figure 2) and water-level fluctuations in observation wells.The conceptual model was designed according to the actual groundwater flow in the basin.It was simulated by three layers in the vertical section and one layered aquifer, where the horizontal and vertical flows between the simulated layers were considered.The model layers, from the ground surface down, are shown in Table 2.The model reflects the presence of semi-permeable overlying and underlying sediments and a water-bearing horizon, tapped by the previously mentioned groundwater sources.As such, the influence of any hydraulic contact with sandy water-bearing strata that overlie the tapped horizon is avoided, because of both insufficient analysis andartificially-created groundwater renewal difficulties, which certainly affect groundwater source performance.The layer contours were geometrized and transformed into the coordinate system of the model based on voluminous borehole data collected in the extended area.Figure 4 shows the three-dimensional hydrogeological model of the sectional view of Šumice.The tapped water-bearing horizon is shown in blue, with production well screens shown in green.As previously mentioned, the overlying and underlying strata of this water-bearing horizon are comprised of semi-permable and impermeable sediments.

HYDRODYNAMIC MODEL
The concept of the hydrodynamic model of the groundwater sources that provide water supply to Kikinda is based on the simulation of three-dimensional groundwater flow.Natural factors, such as the type and characteristics of the represented geological units, the distribution of water-bearing and impermeable units, the seepage characteristics of porous media and the mechanism and regime of groundwater flow, as well as the goal of the task at hand, were of primary importance in selecting the mathematical model concept.A multi-layered model was developed.The hydrodynamic model represented the municipal groundwater source of Kikinda (the source at Šumice and the Jezero Well), but also the extended area including the groundwater sources that service the following factories: MSK (producer of methanol and acetic acid), TM (tile producer), and LŽT -Kikinda (manufacturer of industrial machinery).
The tree-dimensional (3D) finite-difference numerical model for the present study was developed using Modflow [20] with Groundwater Vistas as the graphical user interface [21].The model encompassed an area of 6 km × 6 km and a depth of 300 m.Its orientation was north-south and the discretization 100 m × 100 m in three layers.Grid cells size was refined to 12.5 m × 12.5 m in the area of interest, wherewell density and extraction rates were high (Figure 5).To illustrate the flow field discretization and the derived schematization of the geometric relationships between the lithological elements in the study area, Figure 6 shows the result of schematization in the vertical section.

Hydraulic parameters
The seepage characteristics of the model layers were specified by the values of hydraulic conductivity, storage coefficient, specific storage, and effective porosity.These hydrogeological parameters were assigned as representative values to each cell.The initial values of hydraulic conductivity of the model were those obtained from earlier tests of the production wells.During the course of the present research, no well pumping tests were conducted for the purposes of the hydrodynamic model, given that they would have affected the municipal and industrial water supply.The initial values of the storage coefficient, specific storage, and effective porosity were adopted based on international reports [22][23][24][25], related to the hydrogeological properties of sediments with similar characteristics.

Boundary conditions
The following boundary conditions were specified in the hydrodynamic model of Šumice: head-dependent flux boundary condition (Cauchy or mixed conditions) and boundary of prescribed flux (Neumann conditions).
The head-dependent flux boundary condition (Cauchy or mixed conditions) -the influence of the Galacka Canal, was simulated using this boundary condition.Given the depth of the pumped aquifer at the source and the thick package of semi-permeable and impermeable sediments in its overlying bed, the surface streams have no affect on the groundwater regime of the tapped horizon.This boundary condition was specified in the first layer at locations where both canals exhibited minimum multiannual values, or 74 m.a.s.l.(metres above sea level).The canal widths were specified according to a topographic map, with an elevation 0.5 m less than the minimum canal water level, and the thickness of the bottom canal layer of 0.15 m, with a hydraulic conductivity of 1 × 10 -6 m/s.
Using the same boundary condition (head-dependent flux boundary condition), the effect of the source of recharge/point of drainage, located outside of the model area covered, was simulated.In the Modflow code applied here, it was represented by the "general head boundary".This approach was used to set the registered piezometric levels of the pumped water-bearing horizon in the second layer, given the remote locations where this horizon is recharged.This type of boundary condition is shown in Figure 5.
The prescribed flux boundary (Neumann conditions), or the impact of the wells on the groundwater source, was simulated through a specific flux boundary.In this case, the flux was specified as a function of position and time.Figure 5 shows the positions of the wells, whose operation was simulated in the model with a specified prescribed-flux boundary condition.For the purposes of the hydrodynamic model, well discharges at Šumice were registered every seven days during the period from January 1 to December 31 (Table 2).Figure 7 shows individual well discharges recorded at Šumice during the one-year period.The individual well discharges in the model of Šumice were specified in accordance with both the recorded operating hours and discharge rates.The well screens were specified using as-built dimensions and positions.The discharge rates of the production wells of the other groundwater sources in the area were specified based on available technical reports [26], as constant values throughout the period.These included the MSK groundwater source, operating two production wells with a total capacity of 54 l/s; the TM groundwater source, using a single production well with a discharge rate of 35 l/s; and the LŽT -Kikinda groundwater source, which operated two production wells with a total capacity of 50 l/s.Effective infiltration.In the overall groundwater balance, for the reasons stated in the description of the boundary condition "rivers", the so-called vertical balance had no direct impact on the tapped water-bearing horizon.It was simulated in the first layer by variable-flux Neumann boundary conditions.This value is the sum of infiltrated precipitation and evapotranspiration.The depth-to-groundwater during the entire study period (January 1 to December 31, 2009) was several meters below the ground surface (4-7.3 m).Bearing in mind that the top layer is semi-permeable, the impact of rainfall infiltration on the groundwater regime is small.Ten-percent infiltration was used as the initial value of effective infiltration.Figure 8 shows mean monthly precipitation levels for the period from 1996 to 2006.

Model calibration
The model was calibrated under the conditions of transient flow, with a time step of seven days for the study period (January 1 to December 31, 2009).Figure 9 is a parallel representation of the total capacity of Šumice and the registered piezometric levels in the observation wells.The calibration of the model was deemed completed when a satisfactory match between registered groundwater levels and those obtained by calculations was obtained.Figure 10 shows the distribution of piezometric levels in the tapped water-bearing horizon for the maximum rate of groundwater extraction at Šumice of 117.66 l/s (August 20, 2009).Figure 11 shows the groundwater levels registered at the observation wells and those obtained by calculations in the model calibration process (using the same observation wells).The agreement of the registered and calibrated groundwater levels was rather good.

GROUNDWATER BALANCE
Assessment of the groundwater balance of the area covered by the model revealed that under the conditions that prevailed during maximum groundwater pumping at Šumice, most of the groundwater in the tapped water-bearing horizon came from the south (35.23%).Table 3 shows the inflow into the area covered by the model under the given conditions (August 20, 2009), or, in other words, the summary groundwater balance of municipal (117.66 l/s) and industrial (139.00l/s) sources of water supply.

RADIUS OF INFLUENCE OF THE GROUNDWATER SOURCE AT ŠUMICE
Modpath [12] was mostly used to simulate the conservative advection of dissolved phase contaminants or microbes within the groundwater system over selected periods, ignoring the effect of dispersion.Hydraulic heads from Modflow [20] were used and pore water velocity computed based on hydraulic conductivity and porosity data.The three-dimensional models easily allowed the consideration of layer heterogeneity and partial penetration of the well screens.As previously stated, particle tracking is suitable for different purposes.In this case, the Modpath code was used to assess the radius of influence of the Šumice wells.The distance from which groundwater reached the wells in 200 days was determined.The streamline arrows pointing to the production wells denote the 20-day travel time of a conservative particle.Figure 12 shows that the radius of influence of Šumice (blue contour) is 700-750 m and that of the Jezero Well 350 m.Any paper that does not address new and innovative aspects of the topics of the meeting may be rejected by editors before entering the review process.In order to process the reviewing in time, please submit your manuscript via web interface in camera-ready form and in Microsoft Word format.

CONCLUSION
A groundwater flow model was developed and calibrated for a Quaternary aquifer under unsteady-state conditions.Assessment of the groundwater balance revealed that most of the groundwater flow to the water-bearing horizon of the groundwater source at Šumice came from the south (35.23%).Assuming Darcian unsteady-state three-dimensional flow and based on the major hypotheses of this method (linear form of the sink term in the mass transfer equation and negligible order of magnitude of dispersion effects), the particle tracking method for the determination of the radius of influence of the groundwater source was used and presented.The maximum well discharge capacity of Šumice was found to be 117.66l/s, at which the radius of influence of Šumice was 700-750 m and that of the Jezero Well 350 m.

Figure 1 .
Figure 1.Geographic position of study area

Figure 2 .Figure 3 .
Figure 2. Distribution of wells at the groundwater source of Šumice

Figure 4 .
Figure 4. 3D hydrogeological model of the groundwater source of Šumice, 230 o azimuth section

Figure 8 .
Figure 8. Mean monthly precipitation levels for the period 1996-2006

Figure 9 .
Figure 9.Total capacity of the groundwater source at Šumice and piezometric levels in observation wells (January 1 to December 31, 2009)

Figure 10 .Figure 11 .
Figure 10.Piezometric head distribution in the extended area of the groundwater source of Šumice (August 20, 2009), at maximum capacity

Figure 12 .
Figure 12.Disribution of streamlines around production wells at the groundwater source of Šumice at a discharge rate of 117.66 l/s, which indicates the distance from the well covered in 200 days

Table 1 .
Hydrogeological parameters from previous research

Table 2 .
Model layers and corresponding terrain layers