The Possibility of Managed Aquifer Recharge (MAR) for Normal Functioning of the Public Water-Supply of Zagreb, Croatia

: With its quantities of groundwater, the Zagreb aquifer is an irreplaceable two with three absorption wells upstream of the Mala Mlaka water pumping station (injection of 300 L/s each and 500 L/s each). The most favorable method to recharge artiﬁcially the Zagreb aquifer near the Mala Mlaka pumping station is achieved with an inﬁltration facility using an elevation of 115 m a.s.l. The use of such a facility will enable the smooth operation of the water pumping station and the possibility of increasing the pumping quantities at the Mala Mlaka water pumping station for the future development of the area.


Introduction
The Sava River is the largest tributary of the Danube River. The springing zone is located in Slovenia and has a confluence with the Danube River near Belgrade in Serbia. It flows through Croatia for a large part of its course and forms the border between Croatia and Bosnia and Herzegovina in the eastern part of Croatia. In the wider area of Zagreb, Observing the Zagreb aquifer complex in which permeable and less-permeable layers are present, the vertical component of groundwater flow in some areas is significantly reduced, which is crucial for the protection of the deeper part of the aquifer and determining methods of capturing groundwater. Therefore, the Zagreb aquifer ( Figure 1) is in a hydrogeological sense vertically stratified into two aquifers, separated by a low permeable layer of fine-grained clay sediment [1,2]. Detailed hydrogeological characterization is presented later in this paper.
Previous research on the Zagreb aquifer focused on hydrogeological interpretations and determining the connection between the geological structure and the Quaternary aquifer system of the Sava River in Croatia [3], determining the impact of Zagreb landfills

Mathematical Model of the Zagreb Aquifer
The hydrogeological system of the Zagreb aquifer ( Figure 1) is a complex natural system, which was simplified using a mathematical model for meaningful analysis. The final aim of the mathematical model was to capture the distribution of potentials (piezometric water levels) and to simulate groundwater flow for the purpose of selecting the best MAR solution.
The input parameters for the mathematical model of the Zagreb aquifer were groundwater levels, water levels of the Sava River, precipitation amounts, and pumping quantities for the needs of the public water-supply. Groundwater levels were determined based on data from 392 measuring points (piezometers) in the study area. The water levels of the Sava River in the investigated area were measured at a total of nine gauging stations uniformly distributed from the entrance profile to the model on the western part of the model to the exit from the model in the eastern part of the model. To construct the mathematical model, we used data from four rain gauging stations that covered the model area well.
Prior to the development of the mathematical model of the groundwater flow in the Zagreb aquifer, a conceptual model of the Zagreb aquifer was defined, based on existing research and databases, i.e., geological and hydrogeological characterization, the geometry of the aquifer system, types of the aquifer, porous medium properties, recharge and discharge places of the system, boundary conditions of the model, and initial conditions in the model.

Geological and Hydrogeological Characterization
The wider investigated area is built of Quaternary age deposits in the lowland area along the Sava River in the vicinity of Zagreb. Quaternary age deposits are sedimented in a morphologically irregular area of separated deep basins. In the western part, out of the investigated area, the Samobor aquifer is located with a maximum depth of about 50 m, followed by the Podsusedski threshold where Quaternary deposits are in some places only 3 m thick. After this threshold follows the Zagreb basin, which gradually deepens and in the deepest part in the area ofČrnkovac, it reaches a thickness of about 100 m [1,2], while in the eastern part of Zagreb the aquifers are not as deep. This Zagreb basin is called the Zagreb aquifer. The eastern of the Zagreb aquifer and near the town of Sisak is the new threshold, where the Sisak aquifer is located. The Samobor and Sisak aquifers were not considered in this study, but their citation is important for a comprehensive understanding of the position of the Zagreb aquifer.
The Zagreb aquifer formed during the Middle and Upper Pleistocene when this area was a lake or a swamp area. From the surrounding land area, through the processes of erosion and denudation, the deposits were eroded and washed away by numerous streams into that lake or swamp area. At the beginning of the Holocene, the Sava River breached and then the transport of materials from the Alps began. The water inflow in the hypsometrically lower Zagreb basin led to a drop in water energy and the creation of accumulation conditions with the deposition of a large amount of predominantly wellpermeable gravel up to 60 m thick, above the previously deposited finely grained swamp sediments ( Figure 2). This contributed to the additional vertical separation into two aquifers having different hydraulic characteristics, separated by impermeable or poorly permeable deposits [23,24]. The above-mentioned aquifer geometry was confirmed with a large number of piezometric wells that were drilled in the Zagreb aquifer system. There are also extensive data available on lithology by aquifer depth. The piezometric wells were geologically determined by piezometric well logs data that describe lithology in intervals.
Overall, the hydrogeology complexes were defined using over 1360 piezometric wells and based on 445 selected characteristic piezometers that well reflect the geometry of the aquifer and its lithological relations. It is significant that many of these wells show the entire profile of aquifers, from covering deposits and gravel aquifers to bottom impermeable deposits (clay and marl). In some places, the wells are not evenly distributed because, e.g., the concentration of wells is higher in the locations of existing and planned water pumping stations.
The initial values of the spatial distribution of the hydraulic conductivities for the aquifer complex were obtained entirely based on discrete values because there were no other data. The existing database on hydraulic conductivity was obtained from Croatian Waters (state agency for water management), and lithological data were obtained from individual piezometers [2] or several hydrogeological studies related to determining the hydrogeological parameters of the Zagreb aquifer [23,[25][26][27][28][29]. The database on hydraulic conductivity contains approximately 350 discrete values of hydraulic conductivity, determined as part of water research works in the Zagreb aquifer.
The values of specific yield and porosity for individual aquifers have not been as well researched in the Zagreb aquifer; therefore, the required values were determined based on lithological data from existing piezometric wells and partly obtained from the literature [30][31][32].
The spatial distribution of hydrogeological parameters over the entire model domain was made by interpolating discrete values collected from verified piezometric wells.

Covering Deposits
The surface of the investigated area is covered by a relatively thin humus layer that is 0.2-0.3 m thick, below which is a fine-grained clay layer 0.5-2.0 m thick near the Sava River watercourse; moving away from the river the layer thickness increases to about 10 m and is composed of particles of powder, clay, and fine-grained sand of the Upper Pleistocene to Holocene age. In some places, this poorly permeable layer is missing, and therefore, the aquifer is completely open to the surface.
Those cover deposits are important for the protection of aquifers, but within the investigated area, their spatial distribution does not represent a significant barrier. Hence, it was concluded that the aquifer below them is an open aquifer.
The spatial distribution of hydrogeological parameters of covering deposits was determined so that all cover deposits were assigned the same mean effective value of individual parameters. By averaging the values, the spatial distribution of the parameters of the cover deposits was greatly simplified, which did not significantly affect the groundwater flow simulation itself and had a positively affected the groundwater flow simulation speed. According to available data, the average effective value of hydraulic conductivity of 10 m/day was determined for cover deposits.

Upper Aquifer
The upper aquifer has an uneven thickness that ranges from 5 m on the threshold to a maximum of 60 m in the area of the Petruševec andČrnkovec water pumping stations. Downstream from the investigated area its thickness gradually decreases to 10-20 m.
According to available data, the upper aquifer consists of well-compacted unsorted gravel with a sand content of 20-30% and powder up to 5%. The sediment grain size is larger along the Sava River, whereas toward the east, the grain size decreases and gradually turns into sand. Locally, there are layers and lenses of dusty sand 0.2 to 1.3 m thick.
In the western parts of the aquifer and near the Sava River watercourse, the hydraulic conductivity is over 3000 m/day, while in the central parts of the aquifer, it is around 2000 m/day. In the eastern part of the upper aquifer, gravel is completely absent, and only sand is present. Therefore, the modeled hydraulic conductivity decreases to about 300 m/day [33].

Aquitard
At a depth of about 30 m, a weakly permeable clay layer with an average thickness of about 2.5 m (1-11 m) is spread over the largest part of the aquifer area. This was determined on numerous piezometric wells. However, it is not present in some places; therefore, it is an incomplete layer.
The hydraulic conductivity of this layer was estimated to be K = 1 × 10 −6 m/s. Due to its hydrogeological characteristics, when present, it divides the aquifer system into two parts.

Lower Aquifer
Below the weakly permeable clay layer is a lower aquifer about 20 to 30 m thick. The increased part of fine-grained components, reducing the hydraulic values of the lower aquifer and the groundwater, is also under moderate pressure.
The hydraulic conductivity of the lower aquifer layer was estimated based on research conducted in the wider area ofČrnkovec and Kosnica water-supply sites. In the central part of the lower aquifer, the value of hydraulic conductivity is about 750 m/day; toward the east, it is significantly lowered to less than 100 m/day. This is about 30% less than that of the upper aquifer.

Conceptual Model
The boundaries of the model are defined by the geological structure located between the mountainous area in the north (Medvednica Mt.) and the hilly part in the south, which Water 2021, 13, 1562 7 of 24 stretches from Podsused in the west to Rugvica in the east. A similar model domain has been used in other previous studies of the Zagreb aquifer [23,25].
Based on the previously described Zagreb aquifer complex, the conceptual hydrogeological (hydrostratigraphic) model of the Zagreb aquifer system was defined by the existence of two aquifers separated by a thin incomplete aquitard. The extension, thickness, and hydrogeological properties of these layers within the model domain are inconsistent.
Hydrogeologically, the Zagreb aquifer is an open aquifer, for which the upper limit of saturation is also the free groundwater level. In some places, poorly permeable cover deposits are present. However, they do not significantly affect the groundwater flow in a conceptual sense, but it is still possible that minor anomalies may occur in groundwater flow simulations, especially in the parts of the Zagreb aquifer where the impact of cover deposits is somewhat higher. These parts of the Zagreb aquifer can then be considered semiconfined aquifers.
Mathematical simulations of groundwater flow require defining the hydrogeological parameters of a particular aquifer. In accordance with the conceptual model of the Zagreb aquifer system, we processed the available data necessary for the assessment of the hydrogeological properties of the cover deposits, the upper aquifer, and the lower aquifer. The parameters primarily included the assessment of the values of the vertical and horizontal hydraulic conductivity (K) of individual aquifers, specific yield (Sy) for open aquifers, and the effective porosity (n) of individual aquifers in the Zagreb aquifer. For the needs of the model, the hydraulic conductivity Ky is equal to the hydraulic conductivity Kx, whereas the value of the hydraulic conductivity Kz is based on available data and previous research, taken as 10-times smaller values.
The analysis of equipotential and water levels of the Sava River in the predicted period of groundwater flow simulation in the Zagreb aquifer, i.e., in the period from 2006 to 2010, showed that the aquifers are recharged mostly through the infiltration of water from the Sava Riverbed. The groundwater flow direction is generally from west to east or southeast depending on hydrological conditions. The marginal parts of aquifers are hydraulically impermeable in the north, the main inflow is from the west, and the outflow from the aquifer occurs in the eastern part. Research indicated that there are certain inflows of varying intensity along the southern boundary of the model [34], which was considered during the development of the model.
A certain infiltration of precipitation into the aquifer system should be considered, as well as the underground inflow of water from the upstream area and the inflow from the southern hilly area. Comparing groundwater levels from piezometers and data on water levels from the gauging station on the Sava River, a strong correlation was observed, especially for piezometers closer to the Sava River. As the distance from the Sava River increases, this correlation decreases but is still present.
The above is confirmed by previous calculations of the water balance of the Zagreb aquifer, according to which the contribution of the Sava River to the restoration of groundwater was estimated to be about 73% [35].

Boundary Conditions
The boundary conditions of the water flow model were set in accordance with the conceptual model of the Zagreb aquifer system and additionally defined by the spread of less permeable deposits on the edge parts of the Zagreb aquifer.
The four types of boundary conditions of the model (main inflows and outflows from the model domain) and boundary conditions within the model domain (influence of the Sava River, active water pumping stations, and the impact of the Sava-Odra canal) were defined as shown in Figure 3.  The northern boundary of the model domain (inflows from the direction of Medvednica Mt.) is defined as practically impermeable, although there are numerous, mostly smaller, torrents along it that form in the submountain zone, but are practically all channeled by the infrastructural development of Zagreb. Therefore, they do not have a direct impact on the Zagreb aquifer. If there is any impact on the Zagreb aquifer, then it is mostly smaller and occasional; hence, it would not significantly affect the groundwater flow simulation.
The main groundwater inflow into the Zagreb aquifer is defined in the western part of the model boundary, which is related to inflows from the upstream part of the Samobor aquifer. There are also several smaller inflows along the southern boundary of the model, along which a smaller inflow into the Zagreb aquifer was identified in several places [34]. The main runoff from the system is defined on the southeastern boundary of the model domain, from the Sava-Odra Canal area, across the Odra watercourse, to the Sava River downstream of Zagreb. All boundary conditions of the model domain were simulated by the second type of boundary conditions (Flux-or Neumann-type boundary).
The infiltration of total precipitation in the Zagreb aquifer is defined by the second type of boundary conditions, as a model input over the entire surface of the model domain.

Role of the Sava River within the Boundary Conditions
In the initial simulations of groundwater flow, the Sava River was simulated along its entire length as the first type of boundary condition (head-or Dirichlet-type boundary). Conceptually, this was justified because it is a surface water flow that was assumed to be in strong good hydraulic connection with the Zagreb aquifer; therefore, its effect could be simulated in the form of known piezometric potential. However, shortly after initial groundwater flow simulations, we noticed that its impact is not so simple that it could be simulated only using a Dirichlet-type boundary and that the impact of the Sava River is one of the crucial factors affecting the total groundwater flow.
The Sava River does not represent a watercourse that intersects the entire aquifer system along its entire depth so that the hydraulic connection of the Sava River with the Zagreb aquifer system is not uniform. Therefore, it was ultimately simulated over the entire domain of the Sava River flow model using the third type of boundary condition (transferor Cauchy-type boundary). This type of boundary also represents the dependence of the model on the piezometric potential (water level), but it is also conditioned by the resistance to water flow that occurs during the infiltration of water from the Sava riverbed into the aquifer, and vice versa. This means that the Sava River influences the surrounding aquifer (and vice versa), which is realized through the transition medium, the clogging layer, a certain hydraulic conductivity, and thickness.
Therefore, the conceptual model was corrected and revised, and the Sava River was divided into a total of 34 separate segments. The segments represent a uniform distribution of the entire watercourse and were selected based on previously interpolated data from gauging stations on the Sava River. In the processes of successive simulation and calibration of the mathematical model, transfer in/out relations were determined for each segment. In the end, a satisfactory solution of the Sava River's influence on the water flow in the Zagreb aquifer was obtained.

Water Pumping Station Locations
In the area of the Zagreb aquifer, several water pumping stations for public watersupply have been built for the needs of the cities of Zagreb and Velika Gorica. This large aquifer contains large reserves of groundwater, but Zagreb being located directly on the aquifer and the cover deposits being relatively thin (up to several meters) create difficult conditions for protecting the quality of groundwater in the area. Part of the water pumping station is active, part is spare, and part is, unfortunately, abandoned due to the pollution of significant parts of the Zagreb aquifer.
The Mala Mlaka water pumping station is located on the right bank of the Sava River in the central part of the Zagreb aquifer. It is one of the oldest and largest water pumping stations in Zagreb. Although the idea for its construction started back in 1934, construction began in 1956, and the first wells started operation in 1964. In the first phase, 10 dug wells with a diameter of 6 m and a depth of 13 to 15 m were constructed, which provided about 1.7 m 3 /s of drinking water, which was more than enough for the needs of the city at that time. Subsequently, due to the lowering of groundwater levels, an additional six wells were drilled, which affected deeper parts of the alluvial aquifer and the pumped water is transferred to the dug wells. On average, about 1.3 m 3 /s is produced from the water pumping station.
Data on pumping quantities used in the mathematical model of the Zagreb aquifer refer to the water pumping stations of Mala Mlaka, Petruševac, Velika Gorica, Sašnjak, Žitnjak, and Zaprude, which were active during the modeled period.
Finally, the existing water pumping stations within the Zagreb aquifer were simulated using the fourth type of boundary conditions (singular inputs/outputs within the model). Within the mathematical simulation, this type of boundary condition was separately defined for each well within each water pumping station, considering the position of the well, pumping quantity, the permeability of the surrounding deposits, and the exact position of the filter within the well (multilayered-well definition).

Model Discretization 2.4.1. Spatial Discretization
The horizontal and vertical discretization of the model domain is conditioned by the conceptual model of the Zagreb aquifer and the boundaries of the model. It was constructed made using a FEFLOW finite element network generator, which creates super-elements, i.e., geometric frames of finite elements. In the case of the Zagreb aquifer, three-node triangular elements for 2D model representation and six-node prisms with a triangular base for 3D model representation were selected for discretization. The number of these super elements was determined based on factors that could influence the behavior of the mathematical model in later modeling, such as water pumping station locations, the total size of the model domain and individual super elements, conceptualization of natural conditions, and computer processor speed. Discretization achieves the required heterogeneity, i.e., greater or lesser discretization of individual parts of the model domain. In the case of the Zagreb aquifer, spatial discretization of the model was increased near the Sava River, boundaries of the model, larger lakes, existing active water wells in the Zagreb aquifer, as well as in the wider area of Mala Mlaka and the Sava-Odra Canal ( Figure 4). In the conceptual sense, the Zagreb aquifer is three-dimensionally divided into three separate layers: cover deposits, and upper and lower aquifer.

Temporal Discretization
The temporal discretization of the mathematical model was performed based on the time observation frequencies of three key factors influencing the simulation of groundwater flow in the Zagreb aquifer: Sava River water level, groundwater level at the piezometers, and frequency of pumping volume measured at water pumping stations. The availability, frequency, and uniformity of these measurements in the Zagreb aquifer are not the same everywhere.
The water level of the Sava River is one of the parameters we used in the simulation of the groundwater flow of the Zagreb aquifer. In addition, the available time series of inflows and outflows at the domain boundaries were also considered. The estimates for the inflow and outflow in the mathematical model of groundwater flow are based on the measured groundwater levels near the individual model boundary, i.e., the recorded time series of these measurements.
The groundwater level measuring frequency at the piezometers in the Zagreb aquifer is not uniform. In some piezometers daily measurements are available, in others, there are only two to three measurements per week. Available measurements of pumping quantities at active water pumping stations of the Zagreb aquifer are defined on a monthly basis. Therefore, within the model, these data posed the biggest problem in the temporal discretization of the model. Initially, the time interval was set to one day, but a comparative analysis of the simulation results for a time interval of one day and five days showed no differences. Therefore, a time interval of five days was selected for the final time interval.

Initial Conditions
To solve the nonstationary problem of groundwater flow, it was necessary to define the initial conditions in the model, i.e., the initial distribution of potential (groundwater level) in the Zagreb aquifer for a certain date. The required levels were obtained by searching a database containing time series with measured groundwater levels; in the case of the Zagreb aquifer, these included a piezometric well database and a database of gauging stations on the Sava River.
Measurements of groundwater levels in the Zagreb aquifer began in early 1972 and continue to this day. The series of groundwater levels on several piezometers in the area of the Mala Mlaka pumping station were analyzed. In the period from 1972 to 2019, two characteristic periods with different trends were observed. The first occurred from the beginning of 1972 to 1994, when there was a continuous decrease in groundwater levels in the Zagreb aquifer; after 1994, the groundwater level stabilized, and no further decline was recorded. Continuous decrease in the first period was the result of several different factors. In the initial part of that period, an embankment was built along the Sava River for flood protection, which prevented the flooding of the coastal area but also reduced the infiltration into the aquifer. The second reason could be attributed to the construction of reservoirs and hydropower facilities upstream from Zagreb in neighboring Slovenia. The third reason could be the increasing the amount of pumping rate for public water-supply in Zagreb in that period. All those reasons led to the constant lowering of groundwater levels and some industrial objects, such as thermal power plants, started to have problems with using the water for cooling purposes. Due to that, an artificial river threshold was built in 1993, which caused the stabilizing of the groundwater levels. Looking at the period from 1993 until today, there were no new regulations of the Sava River in the area of Zagreb, and pumping quantities were not significantly increased.
For the period of the mathematical model, a period of five years was chosen (2006-2010) in which there were no large water level variations, i.e., when groundwater levels were low but steady ( Figure 5).
After necessary corrections of the input values and using the Kriging interpolation method, a map of groundwater levels was obtained for the initial time of the simulation of the water flow model in the Zagreb aquifer for 1 January 2006 ( Figure 6).

Model Calibration and Validation
The calibration of a groundwater flow model is a procedure that examines the accuracy of a mathematical model. The modeled and actual groundwater levels at observation points in piezometric wells were compared.
It was necessary to review and process the existing database critically of all available piezometers and their measured water levels, and to remove from the calibration process, all those piezometers for which incomplete or questionable values of measured groundwater levels were observed. Usually, the periods for which the data were most complete and in which there were no extreme water phenomena (extremely low groundwater levels, extremely high groundwater levels, floods, etc.) were selected for calibration, but we wanted to include high and low water conditions, that is, at least one hydrological year.
For calibration purposes, 53 evenly distributed piezometers were selected throughout the Zagreb aquifer, and special attention was paid to the narrower area of the Mala Mlaka water pumping station, where piezometers located next to the wells were selected as calibration points (Figure 7). In the process of model calibration, the biggest challenge was to achieve the best possible match between the model and the actual state, comparing the simulated (modeled) water level values at the control points with the actually measured water levels at these same points. This was largely successfully achieved by adjusting the model parameters during the calibration process: hydraulic conductivity, transfer rate parameter for the Sava River, marginal inflows, the influence of precipitation, and other factors. With each change in an individual parameter, the model was restarted until satisfactory results were achieved. The analysis showed satisfactory results with the coefficient of determination (R 2 ) value of 0.9937 (Figure 8a) and the root mean square error (RMSE) value of 0.3312. same points. This was largely successfully achieved by adjusting the model parameters during the calibration process: hydraulic conductivity, transfer rate parameter for the Sava River, marginal inflows, the influence of precipitation, and other factors. With each change in an individual parameter, the model was restarted until satisfactory results were achieved. The analysis showed satisfactory results with the coefficient of determination (R 2 ) value of 0.9937 (Figure 8a) and the root mean square error (RMSE) value of 0.3312.  Following calibration, the model was validated using field data from 27 December 2010 on the same set of piezometers, and we compared them to the modeled data. The analysis showed satisfactory results with an R 2 value of 0.98 (Figure 8b).

Analysis of Groundwater Levels at the Mala Mlaka Water Pumping Station
The exploitation wells B-1 to B-10 of the Mala Mlaka water pumping station ( Figure  9) were built during the 1960s and 1970s when the groundwater levels of the Zagreb aquifer were up to 5 m higher than today [11]. Then, the filter constructions were defined so that the upper elevations of the filter constructions were 1-3 m below the groundwater level in the dry period. However, the upper elevations of the filter structures nowadays remain in the unsaturated aquifer zone during dry periods as a result of the lowering of the groundwater levels in individual wells, significantly reducing the capacity of these wells. Therefore, at the end of the 1980s, the construction of additional drilled wells of greater depths (B-15, B-16, B-17, B-18, B-19, and B-21) began with the above-mentioned dug wells (Figure 9). From them, water is pumped into dug wells during low groundwater levels because there are pumps in the dug wells that are connected directly to the water-supply system.
Based on the mathematical model of the Zagreb aquifer, each of the exploitation wells of the Mala Mlaka water pumping station was analyzed individually, i.e., the relationship between groundwater levels and the level of filter construction during the research period Following calibration, the model was validated using field data from 27 December 2010 on the same set of piezometers, and we compared them to the modeled data. The analysis showed satisfactory results with an R 2 value of 0.98 (Figure 8b).

Analysis of Groundwater Levels at the Mala Mlaka Water Pumping Station
The exploitation wells B-1 to B-10 of the Mala Mlaka water pumping station ( Figure 9) were built during the 1960s and 1970s when the groundwater levels of the Zagreb aquifer were up to 5 m higher than today [11]. Then, the filter constructions were defined so that the upper elevations of the filter constructions were 1-3 m below the groundwater level in the dry period. However, the upper elevations of the filter structures nowadays remain in the unsaturated aquifer zone during dry periods as a result of the lowering of the groundwater levels in individual wells, significantly reducing the capacity of these wells. Therefore, at the end of the 1980s, the construction of additional drilled wells of greater depths (B-15, B-16, B-17, B-18, B-19, and B-21) began with the above-mentioned dug wells (Figure 9). From them, water is pumped into dug wells during low groundwater levels because there are pumps in the dug wells that are connected directly to the water-supply system. Based on the mathematical model of the Zagreb aquifer, each of the exploitation wells of the Mala Mlaka water pumping station was analyzed individually, i.e., the relationship between groundwater levels and the level of filter construction during the research period 2006-2010. For that period, measured and modeled groundwater levels at these wells were compared.
In the study period, during the dry periods when the groundwater levels were the lowest, several wells of the Mala Mlaka water pumping station experienced problems because their filter construction remained in the unsaturated zone of the aquifer. These are wells B-5, B-6, B-7, and B-8, which are located in the central part of the water pumping station (Figure 9). During the study period, the lowest groundwater level was recorded on 17 September 2007, whereas the maximum level was recorded on 30 December 2010 ( Figure 10).

Results and Analysis of Various MAR Solutions
Based on the mathematical model of the Zagreb aquifer, six different variants of MAR were simulated; one was related with constant potential in the Sava-Odra Canal, three with the recharge from the Sava-Odra Canal with different backwater levels in the infiltration facility, and two with the construction of injection wells upstream of the Mala Mlaka water pumping station. Overall, the variant solutions can be divided into two basic groups: (a) MAR through the bed of the Sava-Odra Canal and (b) construction of injection wells upstream from the Mala Mlaka water pumping station. This paper presents the effects of these MAR variants on the exploitation dug well B-7 of the Mala Mlaka water pumping station, where the negative impact of low groundwater levels of the Zagreb aquifer on the normal functioning of public water-supply was most pronounced.
(a) MAR through the Bed of the Sava-Odra Canal The Sava-Odra Canal was built after the extreme floods of the Sava River near Zagreb in 1964 to protect the city of Zagreb from floods. It was built so that it bypasses the city on the south side and manages the 1000-year high waters of the Sava River of 1510 m 3 /s. The Sava-Odra Canal is rarely used, only when the water levels of the Sava River exceed the level of the dam at the entrance to the Sava-Odra Canal. The projected flow of the Sava River at the time of the canal activation was 1900 m 3 /s, but due to the erosion of the Sava riverbed, the current flow of the Sava at the time of the canal activation was 2350 m 3 /s [36]. Since its construction in 1970 until today, the canal has been in operation only a few times: twice in 1979, in 1980, 1990, and 1998; in the period of 2006-2010, it was used only once in the fall of 2010. To increase the functionality of the canal, it was proposed to lower the dam at the entrance to the canal by ensuring that the canal begins to be active at a flow of the Sava River of 1900 m 3 /s [36].
Several MAR variants in the area of the Mala Mlaka water pumping station were simulated with the mathematical model of the Zagreb aquifer. One of the options is MAR from the Sava-Odra Canal while maintaining a constant potential of 0.5 m, but this requires larger construction works on the overflow structure and the entire length of the canal, as well as transferring part of the Sava River water to the canal.
Due to the need for large amounts of water in the Sava River, the option of building a dam with a height of about 2 m in the Sava-Odra Canal at the height of the Mala Mlaka water pumping station of 113 m a.s.l. was considered ( Figure 13). This dam would ensure the constant presence of water in the canal in the narrower area of the Mala Mlaka water pumping station. The MAR was simulated for water levels in the canal of 115, 114.5, and 114 m a.s.l. which enable a water level in the canal of 1-2 m in the immediate zone of the water pumping station. A water level of less than 1 m in the canal was not considered due to the questionability of the implementation of such recharge. Water levels higher than 115 m a.s.l. were not realistic due to the limitation of the water level of the Sava River at the entrance to the Odra canal during the summer months when the demand for recharge is largest [37].
In the area of the Sava-Odra Canal in the considered section around the Mala Mlaka water pumping station, the bottom of the canal consists of cover deposits of the Zagreb aquifer. According to the data from the lithological database [2], the thickness of these deposits ranges from 1.5 to 3 m, which is important for defining the impact of infiltration through cover deposits and determining the value of the transfer rate parameter in the mathematical model. This part of the Sava-Odra Canal was simulated as the third type of boundary condition, with an assessment of infiltration into the aquifer (Cauchy boundary condition). In the mathematical model, the infiltration volume is calculated over the area, the transfer rate parameter, and the difference between the reference level and the actual groundwater level.
Estimates of the hydraulic conductivity of cover deposits, the clogging layer for the Sava-Odra Canal, are 8.2 × 10 −5 m/s [38], i.e., calculated at about 7 m/day. According to several studies [24,36], the vertical hydraulic conductivity of the Zagreb aquifer deposits is 10 times lower, which is about 0.7 m/day. The thickness of the cover deposits in the pumping zone is a maximum of about 3 m, which produces a transfer in parameter value of about 0.2 day −1 , which was used in the simulations.
The impact of the Sava-Odra Canal under the condition of a constant water level in the canal of 0.5 m did not have as strong an impact on the Mala Mlaka water pumping station as on the wider area, i.e., the relatively wide zone around the canal. The impact on well B-7 of the Mala Mlaka water pumping station is an increase in groundwater levels by about 1.8 m. This well experiences many problems during low groundwater levels because the groundwater level falls below the upper edge of the filter construction, and during the extreme dry periods, the whole filter construction remains dry. Applying the MAR of the Sava-Odra Canal with a constant water level of 0.5 m, the levels would be about 0.5 m above the upper edge of the filter (Figure 14).

(b) Construction of Injection Wells
One of the considered options for artificial recharge is the construction of injection wells upstream of the exploitation wells of the Mala Mlaka water pumping station (Figure 18), which would ensure the self-purification of water before reaching the exploitation wells. Two variants of artificial recharge using a battery of three injection wells with a capacity of 300 L/s and a capacity of 500 L/s were considered. The effect of injection wells on the increase in groundwater levels at the Mala Mlaka water pumping station was examined by simulating the groundwater flow.
The impact of the three injection wells with a capacity of 300 L/s upstream from the Mala Mlaka water pumping station on well B-7 of the Mala Mlaka water pumping station is an increase in groundwater levels of about 1.3 m. Applying the MAR with three injection wells with a capacity of 300 L/s, the groundwater levels would be around the level of the upper edge of the filter structure during the summer dry periods, and with each over pumping, problems would occur due to the drying of the filter structure ( Figure 19). The impact of three injection wells with a capacity of 500 L/s upstream from the Mala Mlaka water pumping station on well B-7 of the Mala Mlaka water pumping station is an increase in groundwater levels of about 2 m. By artificial recharge with three injection wells with a capacity of 500 L/s each, the groundwater levels would be slightly less than 1 m above the level of the upper edge of the filter structure during summer dry periods ( Figure 20).

Discussion and Conclusions
The Mala Mlaka water pumping station is located on the right bank of the Sava River, between the Sava-Odra Canal and the Sava River. The idea of building a water pumping station on this location started back in 1934, construction began in 1956, and the first wells started operation in 1964. Initially, 10 dug wells with a diameter of 6 m and a depth of 13 to 17 m were drilled, which provided 1.7 m 3 /s of water for the water-supply system of the city of Zagreb.
At the time of the construction of the wells, the groundwater level of the Zagreb aquifer was a few meters higher than the upper elevations of the filters, which enabled the smooth operation of the water pumping station. From the 1980s until today, there has been a continuous decline in groundwater levels because of the erosion of the Sava riverbed and the overexploitation of the aquifer. As a result, the groundwater level in some wells of the Mala Mlaka water pumping station is lower than the upper levels of the filter structures, which reduces the capacity of the water pumping station. This initiated the construction of additional and auxiliary drilled wells of greater depth, from which groundwater is pumped and transferred to the dug wells, and further pumped to the water-supply system. Eight auxiliary wells have been constructed with a depth of 31 to 41.3 m, with an elevation of the upper edge of the filter of 96.3 to 101.6 m a.s.l. Without the use of these auxiliary wells, the dug wells of the Mala Mlaka water pumping station would be out of order for most of the year today.
In the immediate vicinity, on the south side of the Mala Mlaka, a water pumping station passes the Sava-Odra Canal, which was built to protect Zagreb from flooding. The high waters of the Sava River enter the Sava-Odra Canal through the overflow structure and lower the water level of the Sava River in the area of Zagreb. From its construction until today, the Sava-Odra Canal has been in operation only a few times. The reason for this is the overflow level, which has changed due to the erosion of the Sava riverbed, so only extreme waters enter the canal. Therefore, the idea of a constant potential (constant water level) in the canal is not a realistic option without major construction modifications of the overflow facility and the entire length of the canal. The amount of water that would be necessary to maintain a constant water level in the canal would significantly affect the water levels of the Sava River during dry periods. Nevertheless, the initial simulation was conducted under the assumption of maintaining a constant potential (constant water level) in the Sava-Odra Canal. This simulation showed satisfactory results regarding the minimum groundwater levels in relation to the upper elevations of the filter structures of the wells. However, it requires extremely large amounts of water that can only be transferred directly from the Sava River. This would create water shortages in the Sava River during the summer dry periods and probably negatively affect the rest of the Zagreb aquifer.
Therefore, additional simulations were performed with the assumption of the construction of a 2 m tall barrier facility in the Sava-Odra Canal at the Mala Mlaka water pumping station at an elevation of 113 m a.s.l. This would allow a maximum water level in the artificial recharge facility of 115 m a.s.l., or 2 m at the barrier site. Since the Sava-Odra Canal has a fall, upstream from the barrier point, the depth of water in the infiltration facility would decrease. This would ensure the constant presence of water in the canal in the narrower area of the Mala Mlaka water pumping station. The water of watercourses, which passes through the canal on the south side of the Sava-Odra Canal and flows into the canal, can also be used to fill the infiltration facility.
To simulate the artificial recharge of the Zagreb aquifer from the Sava-Odra Canal, a mathematical model of the Zagreb aquifer was constructed and calibrated according to the measured values of groundwater levels from selected piezometric boreholes distributed throughout the Zagreb aquifer. In the area of the Mala Mlaka water pumping station, piezometers located next to the wells of the water pumping station were also used for calibration additionally ensure accuracy. Data such as pumping rates, hydrogeological parameters of the Zagreb aquifer, groundwater levels, water levels of the Sava River, and local watercourses, as well as the geometry of the aquifers, were used as input parameters for the mathematical model. All of them were organized into a GIS project that facilitated their search and browsing. The model was calibrated for 2006 because the data for that year were the most complete. After the completion of the calibration process, the simulation was conducted for a period of 5 years (2006-2010), which served as the basis for various models of artificial recharge of the Zagreb aquifer.
Groundwater flow for the period of 2006 to 2010 was simulated for six different variants of artificial recharge (MAR). One assumes a constant potential in the Sava-Odra Canal, three are related to recharge from the Sava-Odra Canal with different backwater levels in the infiltration facility (elevations of 114, 114.5, and 115 m a.s.l.), and two with three absorption wells upstream of the Mala Mlaka water pumping station (injection of 300 L/s each and 500 L/s each).
Only one of the simulations produced insufficient results for the entire simulation period: the artificial recharge of the aquifer with three injection wells with a capacity of 300 L/s each; the other recharge methods provide satisfactory results. This means that in all hydrological conditions through the period 2006-2010, the upper elevations of the filter construction of the wells of the Mala Mlaka water pumping station remain below the groundwater level.
The most favorable method to recharge the Zagreb aquifer artificially is achieved with an infiltration facility using an elevation of 115 m a.s.l. The use of such a facility will enable the smooth operation of the water pumping station and the possibility of increasing the pumping quantities at the Mala Mlaka water pumping station for the future development of the area.