Groundwater Nitrate Pollution Sources Assessment for Contaminated Wellfield

Abstract

Nitrates are one of the most common groundwater contaminants and they come from different sources. The paper presents a study of groundwater quality at Varaždin wellfield in the north part of Croatia. The nitrate concentration at this location has been above the maximum allowed concentration for several decades, which has made the opening of new wellfields costly. Based on the previously developed groundwater flow model, a model that covers the narrow area of the wellfield is developed. The influential zone of the observed wellfield in working conditions is determined. Based on the developed model, the sources of nitrate pollution are located, which can be generally divided into non-point and point sources. Considering the time of groundwater retention in the horizontal flow, it is concluded that the water protection zones are marked following the applicable ordinance. Based on the developed groundwater flow model at the observed wellfield, a simulation of nitrate pollution propagation by advection and dispersion processes is performed. The simulation results point out the location of the poultry landfill as the largest source of nitrate pollution. However, poultry farms, which are located in the influence area of the wellfield, also contribute significantly to the nitrate concentration at the wellfield.

1. Introduction

Quality water supply is one of the necessary preconditions for sustainable development, but in the last few decades, due to more frequent and intensive anthropogenic sources, negative impacts on water supply systems have occurred. On the one hand, there is a constant increase in the need for quality water, and on the other hand, due to the increased emission of harmful substances into the environment, the degradation of groundwater quality is almost inevitable. With the growth of the human population and the increase in the quality of life, the need for drinking water has been also growing, so it was estimated that global water consumption has been increased seven times since the beginning of the 20th century [1]. This leads to the opening of new wellfields or increasing the pumping rate of the existing ones. For sustainable development and further improvement of the quality of life, the protection of groundwater as the main source of drinking water should be one of the main global priorities.

Agricultural production is recognized worldwide as one of the main sources of groundwater pollution; primarily because nitrates are an important component of mineral fertilizers. Several studies have shown that the application of fertilizers and pesticides in agriculture significantly increases nitrate concentrations in surface and groundwater [2,3,4,5]. In order to increase agricultural areas and to improve conditions for food production, different amelioration measures are applied, which significantly contribute to the degradation of groundwater quality [6,7,8]. Namely, every lowland catchment presents a sensitive eco-system with a slow flow and a high level of groundwater [9,10]. Such hydrogeological conditions enable pollutants to penetrate easily from surface to subsoil, where they can stay for a long period of time. This violates the natural balance between water and nutrients [11] and increases nitrate concentrations in surface and groundwater [12]. Additionally, the results of numerous studies show that the maximum allowable concentration (MAC) is often exceeded [13,14,15].

In order to protect water resources and to mitigate negative consequences of nitrate pollution resulting from agricultural practices, the European Union has adopted certain measures in the Nitrate Directive document [16]. As a full EU member, the Republic of Croatia is obligated to follow its directions, such as those regarding the protection of its rich water resources [17,18]. In Croatia and the EU, the maximum allowable nitrate concentration (MAC NO₃⁻) in water for human consumption is 50 mg/L NO₃⁻ [19,20].

Numerical models of groundwater flow are now common techniques used for water resource management and they present a simplified version of the real system and the processes that take place within it [21,22]. The models define groundwater levels and groundwater retention time. Numerical models of pollution propagation, which are an upgrade to the numerical model of groundwater flow, are also common techniques today for the purpose of estimating the pollution propagation and searching for their sources. Lee et al. [23] used groundwater modeling to identify the source and transport path of pollution that was caused by inorganic matter and artificial sweeteners in the groundwater of an agricultural area in Korea. Wang et al. [24] conducted a comprehensive evaluation model of ammonia pollution trends in a groundwater source area along a river in residential areas. The simulation results show that the ammonia pollution plume over time mainly affecting the groundwater quality. These methods are also used as tools to improve the accuracy of recognizing zones susceptible to water pollution in agriculturally rich valley of a lowland river [25]. Additionally, a recent research by numerical methods shows how land use and farming practices may affect groundwater quality [26].

The obtained data from numerical modeling are also used to assess the risk of contamination of drinking water sources, to determine the priorities of remediation works, and to optimally define water protection zones in order to preserve or even improve the quality of groundwater at wellfields. Additionally, sources of pollution can be successfully determined using numerical methods. Recent research shows that cesspits and sewer networks can be major sources of nitrate pollution, and the process of reducing concentrations naturally to achieve water potability again takes several decades [27].

The observed location used to be the largest wellfield in the Varaždin‘s water supply network. Due to the high nitrate concentrations, the wellfield has been inactive since 2004. It is believed such a high nitrate concentration at the Varaždin wellfield has been caused by the application of mineral fertilizers in agriculture [28]. However, there are other wellfields in the vicinity of the observed area where agricultural land is also present, but groundwater nitrate concentrations are significantly lower [29]. If the application of mineral fertilizers had been the dominant pollution source at the Varaždin wellfield, the nitrate concentration would have been similar at the other wellfields in this area, too. The research into the causes of such a high nitrate concentration at this particular location has been conducted, and its results imply that mineral fertilizer application is not the only dominant cause for such a condition. Therefore, this paper aims to:

• Based on the previously developed groundwater flow model given by other authors, develop a new model that covers the narrow area of the wellfield and define dominant nitrate pollution sources;
• Develop a numerical simulation of nitrate pollution propagation in groundwater at the observed wellfield;
• Assess contributions from located point pollution sources;
• Point out a necessity that every wellfield with a significant increase in a pollutant concentration should be analyzed separately.

2. Materials and Methods

2.1. The Study Area

The study area is located in the north of the Republic of Croatia while the micro-location covers the wider area of the Varaždin wellfield (Figure 1), located near the town of Varaždin. The wellfield with two other wellfields (Bartolovec and Vinokovščak) is located around the town and the related infrastructure constitutes the regional water supply network of Varaždin.

Groundwater is pumped from the same alluvial aquifer system built of Quaternary sand and gravel and it was found to be of good quality at all the wellfields, except the nitrate concentration at Varaždin wellfield (

Table 1. Arithmetic means of certain groundwater quality parameters at the observed wellfields in the period 1994–2003.

Parameter	Unit	Varaždin	Bartolovec	Vinokovščak	MAC *
Nitrate	mg/L NO3−	79.41	23.96	18.23	50
Nitrites	mg/L NO2−	0.00	0.00	0.00	0.5
pH value		7.45	7.40	7.41	6.5–9.5
KMnO4 consumption	mg/L O2	0.87	0.71	0.79	5
Electro-conductivity	mS/cm	585	560	550	2500
Chlorides	mg/L Cl−	17.64	23.74	14.51	250
Amonium ion	mg/L NH4+	0.02	0.02	0.03	0.5

* MAC–Maximum Allowable Concentration.

). The value of this groundwater quality parameter was higher than at the Bartolovec and Vinokovščak wellfields [29]. In the Varaždin regional water supply network, the nitrate concentration was reduced by increasing the pumping rates at the wellfields Bartolovec and Vinokovščak where the values were below MAC NO₃⁻ [31]. Due to the constantly high level of nitrate concentration, the Varaždin wellfield has been inactive since 2004.

2.2. Geological and Hydrogeological Conditions

The geological structures consist of coarse-grained sediments of gravel and sand (middle and upper Pliocene age), while the underlying sediments are composed of silt and clay (upper Tertiary age) [32].

The aquifer is parallel with the Drava River and its thickness increases from the west to the east. At the far west end, the thickness is 5 m, while in the town of Varaždin it is about 75 m. In the eastern part the thickness of the aquifer increases to its maximum value of 148 m and then it decreases gradually. The aquifer is composed of two permeable gravel and sand layers, divided mostly by a semipermeable silty and clay aquitard [33]. The upper layer of the aquifer is composed of coarse-grained gravel and sand, while the lower part is composed mostly of finer-grained gravel with sand. The groundwater flow is in the west-east direction and it is parallel to the Drava River. The quality of groundwater is significantly higher in the lower layer of the aquifer at the Bartolovec wellfield [29].

The hydraulic conductivity values change slightly in the plan from 100 to 300 m/day, so it is higher in the west and lower in the east what is the consequence of the size drop of the granulometric fractions in the west-east direction. At the location of the Varaždin wellfield, the hydraulic conductivity values were measured to be in the range 185–241 m/day for the upper and 92.6 m/day for the lower layer [30].

2.3. Experimental Data

The groundwater quality at Varaždin wellfield was controlled at eight wells (B3 to B10) and five piezometers (P23, P25, P26, P29, and P30) from the upper aquifer layer by regular sampling in the period from 1994 to 2003. The periodicity of the sampling was at least once a week for the active wells and at least once a month for the piezometers according to the valid regulations [17]. The activity of wells during the observed period had been changing, which means that different combination of wells was active in a particular period. Since the water was drawn from the upper layer of the same part of the aquifer at all the sampling locations, the results of the analysis relate to it.

2.4. Groundwater Flow Model

In the case of a steady groundwater flow in a heterogeneous aquifer in which the main directions of anisotropy coincide with the coordinate axes it can be expressed by the following expression [34]:

$$\frac{\partial }{\partial x}\left(K_{xx}\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yy}\frac{\partial h}{\partial y}\right)+\frac{\partial }{\partial z}\left(K_{zz}\frac{\partial h}{\partial z}\right)+q_s=0$$

(1)

where

K_ii 

- hydraulic conductivity tensor (L/T);

h 

- groundwater level (L);

q_s 

- volumetric flow representing the source/sink of water in case of pumping/extraction water in/out of the system (L³/T);

The solution of Equation (1) was carried out numerically by the finite difference method and for this purpose, the MODFLOW code was used [35] as part of the software tool GMS 10.3 (Groundwater Modeling System). The simulation was performed for steady flow conditions.

2.4.1. Model Area

The domain of the model is determined in such a way that a sufficient distance in all directions from each well, where the effect on groundwater flow under working conditions of the wellfield is negligible is taken. Namely, some research shows that the well influence radius in unconfined aquifers with a hydraulic conductivity value of 200 m/day can be expected up to 1500 m [36], while the Kusakin [37] and Siechardt [38] methods show the well influence radius for the observed wellfield below 1000 m. As recent research suggests careful use of the well influence radius [39], a value of 3000 m as the minimum distance from the well to the edges of the domain is used. In the simulation process, it is found that the change of groundwater level at the domain boundaries due to a different amount of pumping rates at the wellfield is negligible, which justifies described approach.

The horizontal discretization is done by an orthogonal grid using the finite difference method. The size of the spatial elements is 50 m × 50 m, to be reduced to 25 m × 25 m in the area of the wellfield, due to a more accurate display of the flow in the areas of larger gradients. The vertical discretization is made using four layers representing the aquifer top deposits, then the upper aquifer layer, the low-permeability layer, and the lower aquifer layer (Figure 2). The thicknesses of aquifer layers were taken from the results of previous researches [30,32].

2.4.2. Hydrogeological Model Parameters

Values of hydraulic conductivity were obtained in the calibration process of groundwater level given in the previous study, where values of 213 m/day for the upper and 92 m/day for the lower aquifer layer at the observed area were determined [30]. In the model area, horizontal isotropy was assumed, while vertical anisotropy was expressed by the ratio K_h/K_v = 1/10. The values of the horizontal and vertical hydraulic conductivity of the aquifer top deposits (silt and sand) are 0.1 m/day and 0.01 m/day, and the values of the low-permeability layer of clay, which divides the aquifer into two parts are 1 × 10⁻³ m/day and 1 × 10⁻⁴ m/day, which corresponds to the soil type value range in previous researches [40]. The values of effective porosity are estimated based on the lithological composition of the aquifer, which for both aquifer layers are 23% [41], while for the top deposit and low permeability layer are estimated to be 38% [42].

2.4.3. Groundwater Flow Model Setup

At the Varaždin wellfield, the pumping rate was systematically reduced from 1998 to 2004, when the wellfield was set to inactive mode. To determine the impact of the wellfield in full working conditions, the value of the pumping rate before the pumping reduction (the year 1997) was simulated. Eight wells from the upper aquifer layer were active (B3–B10) with a total annual pumping rate of 9,183,429.20 m³. To simulate the pumping rate at the wells, the daily average of 3145.01 m³/day was calculated. This value was used for each well to determine the impact area of the wellfield. Within the GMS software, theNeumann boundary condition–well is used in this purpose.

Groundwater precipitation recharge is imposed on the surface cell layer, which represents the boundary with a constant inflow. To determine the vertical balance, the precipitation regime at the hydrometeorological station of Varaždin (Figure 1) in 1997 was considered. In that year, the total rainfall was 880 mm. Based on previous research in the wellfield area, the value of effective infiltration of 35% of total precipitation (308 mm) was adopted [43]. Within the GMS software, the Neumann boundary condition–precipitation is used in this purpose.

Dirichlet boundary condition–groundwater head is set at the model boundaries while the exact values of the boundary heads are obtained through the calibration process. This boundary condition allows the determination of the groundwater level and the groundwater flow is simulated at a speed sufficient to satisfy a given gradient. Calibration is completed using the measured values at the control points shown in Figure 1 [30].

2.5. Nitrate Pollution Propagation

Due to their high solubility, nitrates tend to be carried by water. Gardner [44] concluded that two main processes are involved in nitrogen transport: advection and dispersion. The mathematical description of the transport of pollution in groundwater, without chemical reactions and decomposition of pollutants, can be expressed by the advection-dispersion equation [45]:

$$\frac{\partial \left(nC\right)}{\partial t}=\frac{\partial }{\partial x_i}\left(n{D}_{ij}\frac{\partial C}{\partial x_j}\right)-\frac{\partial }{\partial x_i}\left(n{v}_{i}C\right)+q_s{C}_s$$

(2)

where:

C 

- pollution concentration (M/L³);

n 

- porosity;

t 

- time (T);

x_i 

- distance along the corresponding Cartesian coordinate axis (L);

Dij 

- hydrodynamic dispersion coefficient tensor (L²/T);

v_i 

- groundwater flow rate (L/T);

q_s 

- flow per unit volume of aquifer representing the source (positive) or sink (negative) of fluid into the aquifer (1/T);

C_s 

- source or sink pollution concentration (M/L³).

The solution of Equation (2) is carried out numerically by the finite difference method. For this purpose, the MT3DMS code [45] was used as part of the software tool GMS 10.3 (Groundwater Modeling System) and it represents an upgrade to the groundwater flow model using the MODFLOW finite difference method.

The simulation was performed as advection-dispersion transport of nitrates. Although the denitrification process can take place in areas with higher proportions of fine-grained sediment and organic matter in the aquifer, previous research shows that these processes are very weak or even insignificant and that nitrates are transported over long distances by groundwater [46]. Having in mind that the results and the fact that the model domain in this study is relatively small, the approach of nitrate transport to the processes of advection and dispersion is justified.

2.5.1. Hydrodynamic Dispersion

The coefficient of hydrodynamic dispersion is defined as the sum of effective molecular diffusion and mechanical dispersion. The role of molecular diffusion in groundwater pollution transport is very small and can be neglected in modeling pollution transport [45]. Due to this mentioned reason, this process has not been considered in this study.

Dispersivity is a transport parameter and it often changes with the distance over which pollution travel. A good first approximation for the longitudinal dispersivity is to use a value of 10% of the transport distance [47], while Mallants et al. [48] provide the best estimate for the longitudinal dispersivity based on a few experimental studies and empirical formulas given by different authors [49,50,51]. Having in mind various dispersivity relations given by different authors, the resulting empirical formula accounts for all the empirical investigations given by the three mentioned authors. In this way, the new empirical formula can be incorporated into the numerical model providing authors the certainty that relevant scientific investigations are incorporated into the model averaging the experimental curves and data. The curves created by Neuman [49] and Xu and Eckstein [51] have two common points, so they are taken as two firm points in the construction of the new curve. They represent experimental laboratory and field dispersivity values in the range up to 250 m distance. Dispersivity values for distances greater than 250 m are obtained using numerical values given by Gelhar [50] in a way that the average point was calculated based on the numerical values grouped inside the range of 280 m to 1000 m distance range. The fourth point was calculated as the average point of Gelhar numerical values grouped in a distance range 2000–10,000 m. The linear regression through the four mentioned points produces the following linear equation:

$${\alpha }_L=0.0267x$$

(3)

The value of longitudinal dispersivity was determined based on the distance of the largest point source of pollution to the geographical center of the wellfield using Equation (3). To quantify the effect of three-dimensional dispersivity is necessary to determine the values of transverse horizontal (α_TH) and transverse vertical dispersivity (α_TV). These values were obtained in the model calibration process in a way that the deviation between the measured and simulated values of nitrate concentrations at the control points was minimum. Based on the presented facts, the value of longitudinal dispersivity of 30 m, transverse horizontal 12 m, and transverse vertical 1 m was adopted.

2.5.2. Pollution Propagation Boundary Conditions

The inflow of nitrate into the aquifer is achieved by infiltration through the unsaturated zone. The amount of nitrate in the leachate depends primarily on anthropogenic sources on the surface. In the area of the Varaždin wellfield, there are several anthropogenic sources of pollution that affect the concentration of nitrates in groundwater, which can generally be divided into a point and non-point sources. Non-point sources include agricultural areas located north, west, and south from the wellfield, the urban area of the town of Varaždin located northwest and west, and the natural background level. Point sources primarily include several poultry farms, which are located in the western area, and a poultry landfill on the southwest side of the wellfield.

Using statistical and geostatistical methods in the previous study [52], the nitrate concentration of non-point pollution sources in the Varaždin area was estimated to be 35.92 mg/L NO₃⁻ and this value is set on all boundaries. Having in mind the groundwater flow direction, only the values on the western boundary are of practical importance for the model results in the Varaždin wellfield area.

3. Results and Discussion

The conducted scientific research is based on the findings published in the scientific and professional literature on nitrate pollution of groundwater, as well as in the field of numerical modeling and statistical methods. Having in mind that wellfields are particularly vulnerable due to the impact on groundwater flow, today more attention is paid to the protection of their influencal zones, which is crucial for good quality water [53,54]. The U.S. Environmental Protection Agency defines this zone as as the surface and underground area surrounding a well or wellfield through which pollution is likely to move towards the well or wellfield. The scientific contribution of this research is based on the fact that using numerical models of groundwater flow and pollution propagation in interaction with statistical methods, with knowledge or good approximation of input parameters, point sources of nitrate pollution in the influential wellfield zone can be located and quantified.

The first step in nitrate pollution sources assessment is to obtain the spatial distribution of nitrate concentration at the Varaždin wellfield. For this purpose, nitrate concentration annual averages at wells were observed. A 10-year period from 1994 to 2003 was set when the wellfield was active. During the observed period, the nitrate concentration was higher at some wells compared to other ones. The lowest nitrate concentration was recorded at well B3 and the highest at well B10, while the values at other wells were within that interval (Figure 3). Good consistency of the annual average curves was also observed, so as the concentration at one well increases, it increases at the other wells and vice versa. This confirms the value of the coefficient of correlation of annual average nitrate concentrations at wells B3 and B10 which is calculated to be 0.86.

The analysis of the time series determines the existence of regular ratios of nitrate concentration in wells. For this reason, it is justified to calculate the arithmetic means of the nitrate concentration ($\overline{c}$) during the observed period at each well and piezometer (

Table 2. Arithmetic means of nitrate concentration at control points and deviation values obtained by nitrate pollution propagation simulation.

Well/Piezometer	$\overline{c}$ (mg/L NO3−)	Δ	Δ2
B3	74.14	5.18	26.83
B4	75.35	0.98	0.96
B5	76.70	−3.08	9.49
B6	80.04	−0.08	0.01
B7	76.71	−0.60	0.36
B8	81.42	5.34	28.52
B9	85.92	−0.67	0.45
B10	87.80	0.39	0.15
P23	65.03	3.18	10.11
P25	71.67	0.40	0.16
P26	77.40	−0.18	0.03
P29	68.62	0.44	0.19
P30	87.00	0.47	0.22
∑Δ2			77.48
${\overline{\Delta }}^2$			5.96

). The obtained nitrate concentration values will be used as control values when calibrating the nitrate pollution propagation simulation.

By the developed conceptual model, which is defined by aquifer geometry, hydrogeological parameters, boundary conditions, and calibrated based on previous research, a simulation of groundwater flow in the narrow area of the Varaždin wellfield, is obtained.

In Figure 4, a diagram shows the relation of the groundwater heads at the boundary and simulated values of groundwater levels in control points. Model verification is carried out as presented in [55], so the residual mean error is calculated to be −0.023 m, absolute residual mean error is 0.153 m, root mean square error is 0.172 m, normalized root mean square error is 2.15%, while the value of the coefficient of determination is found to be 0.99.

To determine the inflow zone of the wellfield, a simulation of particle tracking is performed. The MODPATH program code [56] is used for this purpose and it is often used in groundwater resource management today [23,25]. The program uses a semi-analytical particle tracking scheme that allows obtaining an analytical expression of the particle path within each cell by the finite difference method. The particle path is calculated by tracking a particle from one cell to another, assuming that all particles move at an effective water flow rate until they reach the model boundary or satisfy another completion criterion. To determine the inflow area of the wellfield, backward particle tracking where the particles are tracked from the endpoint, in this case from the wells in the reverse direction towards the model boundaries is simulated. The limit of the simulation was set at 4000 days.

The simulation of particle tracking shows that the inflow zone of the wellfield extends slightly southwest and approximately 1 kilometer to the north and south (Figure 5). Additionally, the simulation results show that the water protection zones were marked following the valid national regulation [57]. The estimated retention time of groundwater in the horizontal flow from the boundaries of the first water protection zone to the closest well is estimated to be about 200 days, and from the second water protection zone to the well about 400 days. On the west side of the inflow area, there are mostly agricultural areas, which can be observed as non-point pollution sources. In addition, three poultry farms are located in the inflow zone, one on the south side (PF1), two on the northwest side (PF2 and PF3), and a poultry landfill on the southwest side, which are classified as point sources of pollution (Figure 5). Previous research shows that leachate from poultry farms contains nitrate concentration values of up to 180 mg/L NO₃⁻ [58].

The performed simulation of nitrate pollution propagation is steady, which means that the nitrate inflow is constant in the observed period. The values of nitrate concentrations in the area of the poultry landfill and poultry farms, located in the inflow area of the wellfield, were determined in the model calibration process in such a way that the sum of square residual values between simulated and measured values at control points (B3–B10 and P23–P30) is the smallest. The values of the residuals (Δ) and the squared residuals (Δ²) of the simulated and measured values at the control points are shown in Table 2, while the residual mean error is 0.905 mg NO₃⁻/L, absolute residual mean error is 1.615 NO₃⁻/L, root mean square error is 2.441 NO₃⁻/L, and normalized root mean square error is found to be 13.28%.

A comparison of the measured and simulated values of nitrate concentration at the control points is shown in Figure 6. The correlation coefficient concerning the direction of perfect fit is 0.94, which represents a strong correlation between the two data sets. A t-test was performed and the obtained value of t = 9.14 is higher than the critical value t_α = 2.20 at the 95% significance level, so the model can be considered reliable.

Considering advection–dispersion processes and initial conditions, nitrate concentration values in groundwater in the area of point sources of pollution are estimated and the results are shown in

Table 3. Estimated values of nitrate concentration in the area of point sources of pollution.

Point Source of Pollution	Estimated Value of Nitrate Concentration (mg/L NO3−)
Poultry landfill	114
Poultry farm (PF1)	80
Poultry farm (PF2)	86
Poultry farm (PF3)	87

.

In the upstream western area from the location of the poultry landfill, there is a PDS7 piezometer (Figure 7) from which groundwater in the upper aquifer layer was sampled in the period 2004–2006. During this period, the nitrate concentration at this measurement site did not exceed the value of 41 mg NO₃⁻/L [30]. These results show that in the west upstream area from the landfill the nitrate concentration is lower compared to the east downstream area.

The simulation results indicate that the located point sources of pollution have the greatest influence on the nitrate concentration in the area of the Varaždin wellfield (Figure 7). The estimated initial value in the area of poultry landfills is 114 mg/L NO₃⁻, while in the area of poultry farms it ranges from 80 to 87 mg/L NO₃⁻.

Considering that there are several poultry farms (PF4–PF9) in the wider northwestern area of the Varaždin wellfield, a simulation of nitrate transport from these locations is also conducted. The estimated nitrate concentration value from the poultry farm PF3, as the highest estimated value is used at these locations. The simulation results show that the pollution from these farm locations spreads north and they do not affect the nitrate concentration in the area of the Varaždin wellfield (Figure 7).

The performed analyses indicate that point sources of pollution in the vicinity of Varaždin wellfield are the dominant sources of nitrate pollution, with the focus on the southwestern area as indicated in [59]. Additional simulation to estimate the time required for the pollution to be carried by the advection and dispersion processes from the located point sources of pollution to beyond the boundary of the wellfield was performed. Nitrate input from the area of point sources of pollution was simulated only as an impulse condition, while in the rest of the simulation process there is no nitrate input.

Figure 8 shows the state of the simulation after 100 days. The time required to lose the impact of the poultry landfill when that source would not be active is about 1000 days (Figure 9), while the time required to lose the impact of the furthest poultry farm in the northwestern part of the wellfield is estimated to 3200 days (Figure 10).

The results in this paper justify the presented approach and the methodology can be applied to the locations of other wellfields with high pollution concentration, as an effective method in locating and quantifying pollution sources. However, each wellfield should be approached individually in order to properly setup the research concept and determine the input parameters to obtain reliable results.

4. Conclusions

To determine the impact of the inflow zone of the Varaždin wellfield, a numerical model of groundwater flow was developed and a significant impact of pumping on groundwater flow was observed. In the wellfield influence zone, agricultural areas as non-point sources of nitrate pollution and poultry farms and a poultry landfill as point sources of pollution were located.

Considering the time of groundwater retention in the horizontal flow, it is concluded that the water protection zones are determined according to the applicable national regulations.

The conducted simulations show that it is crucial to locate all point sources of pollution in the vicinity of the wellfield, to obtain reliable quantification results because the absence of only one source significantly reduces the accuracy of the simulation. The results indicate that the largest point source of pollution is the poultry landfill, and a significant contribution is also made by poultry farms located in the inflow area. However, to finally confirm these results, additional extensive field measurements are still needed in terms of groundwater sampling in the immediate upstream and downstream areas of point sources of pollution.

In order to improve the quality of groundwater at the Varaždin wellfield, it is necessary to remove the poultry landfill and regulate the operation of poultry farms, which are located in the inflow zone, to release as little nitrate pollution as possible. If possible, a solution for relocating these facilities should be found in the long term period. Namely, although the protection zones are marked following the valid regulations, the results in this paper show that to improve groundwater quality at the location of Varaždin wellfield it is necessary to further protect the inflow zone and reduce the number of point anthropogenic pollution sources.