Methodology for Determining the Die-Off Coefficient of Enterococci in the Conditions of Transport through the Karst Aquifer — Case Study : Bokanjac – Poličnik Catchment

This paper presents the methodology for determining the die-off coefficient of faecal indicator bacteria (enterococci) when transported in a karst environment. The main problem in exploring karst environments, which this methodology strives to cope with, is lack of field measurements, poor data on karst rock formation, fractures and channels within it, and groundwater level dynamics. The analysed karst catchment (Bokanjac–Poličnik) is situated in the hinterland of the city of Zadar (Republic of Croatia) and covers an area of 235.07 km2. In the water supply wells within the analysed catchment, a frequent occurrence of enterococci was observed. The proposed methodology consists of two basic steps. Preliminary analyses as the initial step were used in the accumulation of certain assumptions related to the detection of increased concentrations of enterococci as well as in determination of the potential source of pollution. In the second step, the analytical model was constructed with the aim of resolving processes of sorption and die-off and determining the dominant factor in the process of natural removal of enterococci when transported in karst environment. Within the model, two parts of the pollutant transport are integrated: vertical percolation and horizontal seepage flow and transport. The mean value of the total die-off coefficient by transport through the unsaturated zone in the analysed case is ktot = 8.25. Within the saturated zone the total die-off coefficient ktot is within the limits of 0.1 and 0.5.


Introduction
Presence of enterococci in the pumping sites for drinking water supply of Bokanjac-Poličnik karst aquifer is an indicator of faecal pollution which is very dangerous for the environment and people living in that area.The Bokanjac-Poličnik aquifer is burdened by pollution caused by seepage of polluted water from the surface (septic tanks).The die-off coefficient of faecal indicator bacteria, such as enterococci (ENT), in typical conditions of groundwater flow through the karst aquifer (fracture porosity) is unknown.The objective of this study is to define appropriate methodology for determining the die-off coefficient of faecal indicator bacteria, ENT, when transported in a karst environment.Besides that, additional goals are to determine self-purification capacity (sorption and die-off expressed with associated coefficients) in case of incident events or continuous pollution processes and to determine the dominant process in removal of ENT in the observed karst environment.The proposed In a karst environment, groundwater flow is mostly conducted through fractures and karst channels where quite high flow velocities can be achieved.Although the matrix can be very porous due to small primary permeability, it has no contribution to the flow to a large extent, which is why diffuse flow happening inside of it, is often neglected.On the other hand, the retention capacity of the matrix is quite big and inside its huge amounts of water can be stored; that water is gradually seeped towards springs, fractures and karst channels.Therefore, karst rock and aquifer are characterised by triple porosity (the matrix, fractures, karst channels), laminarly-diffuse flow inside the matrix and turbulent flow inside the karst channels and fractures.
Different types of mathematical and numerical models could be applied to deal with hydrological and ecological problems of karst aquifers with dual and/or triple porosity depending on available data and identified problems.Distributed parameter models are often applied to systems of karst aquifers [6][7][8][9].Those types of models are usually applied to places where enough data about the features of aquifer is gathered.Those features include permeable and non-permeable boundaries, zones of inter-layer seepage, positions of faults, springs and sinks.Distributed parameter models require division of the aquifer into small representative three-dimensional volumes such as finite difference cells, or finite element cells.Basic hydraulic features of aquifer (hydraulic permeability, transmissivity, storage, effective and total porosity) are spaced and can be distinguished and changed for every model cell or element.Various numerical strategies have been used in karst aquifers analysis over the years.The simplest approach is the one which assumes that karst aquifer conducts oneself equivalently towards the medium of intergranular porosity in the conditions of flow and transport.The approach is called single continuum porous equivalent (SCPE) model.SCPE is based on the acquisition of the continuum scheme and flow simulation by the potential flow equation, and the model itself simulates laminar flow.Relevant research in [10][11][12][13][14][15][16] are just a few examples of successful implementations of the SCPE approach when it comes to researches of water resource.The approach is based on the fact that higher values of hydraulic conductivity are given to the model cells for which it is known or assumed that go through the fractures compared to the surrounding cells which therefore represent a weakly permeable matrix.Computer codes mostly used are different variants of the MODFLOW software.
The research which dealt with controlled discharge of sewage water from septic tanks showed the importance of enough thickness of the unsaturated zone.It also showed that the concentration of the faecal indicator bacteria is largely decreased with the flow heading towards the depth [17,18].However, the research also indicated that the transport of the bacteria through the soil intensified due to more pronounced precipitation, and the increased concentration of bacteria on the pumping sites coincided with the periods of intensive precipitation [19,20].In case that bacteria reach the saturated zone relatively quickly and in high concentration, it will represent danger for the ecosystem in which they are found.Several studies [21][22][23][24][25] have shown that the microorganisms and colloids can easily be transported by preferential flows.Firouzi et al. [26] analysed transport and bacteria removal in a limestone soil under the saturated conditions.He noticed that the breach of bacteria mass is enlarged with the increase of velocity of groundwater flow.Potential presence of pathogens in groundwater could be discovered by monitoring the presence of faecal indicator bacteria such as ENT.A general model of bacteria transport consists of advection, dispersion and the process of attachment (sorption), detachment (desorption), filtration, blocking, die-off in contact with the porous medium and die-off in a liquid phase [27].A high rate of attachment and removal of bacteria on limestone grounds could be ascribed primarily to the fact that calcium carbonate is mostly positively charged when in neutral pH level surroundings which makes it convenient for the attachment of negatively charged bacteria.According to [27], blockage and detachment do not have a significant effect on bacteria transport, and the most plausible reason is the die-off that is happening on the surface of particles, that again enables the particle surface to be free for attachment.The bacterial die-off rate in a liquid phase is under the influence of factors such as temperature, light, pH level, amount of poisoned substances, dissolved oxygen and protozoa.The die-off coefficients, k, of the faecal indicator bacteria are defined mainly for the conditions of flow through the aquifer with the intergranular porosity, and with somewhat different values depending on all the above-mentioned factors [28][29][30].One of the most important factors that affect die-off of the bacteria is the temperature.John and Rose [31] compared approximately ten studies which dealt with the investigation of the die-off of the Escherichia coli (EC) and ENT in a temperature range from 3 to 37 °C.In a total sample of 35 values, the die-off coefficient of the EC varies from 0.01 to 1.5 day −1 .Higher levels of die-off (above 0.5 day −1 ) is noticed with temperatures above 20 °C.According to [27], despite heterogeneous conditions during different experiments, from the gathered data, it can be concluded that the value of the die-off coefficient increases with the temperature rise.The average die-off coefficient value with the temperature of 10 °C is 0.15 day −1 that is 0.5 day −1 with the temperature of 20 °C.Impact of the soil to the die-off rate depends primarily on the local conditions on the catchment, and the effect can be positive or negative which depends on the presence of nutrients, toxic substance, oxygen, etc. Heavy metals have toxic impact on bacteria and with its presence in ground water and soil, die-off level is increased [32,33].The level of die-off of the faecal indicator bacteria in oxygen-saturated water is significantly lower in relation to water with low oxygen saturation [32].
Based on the above mentioned, one can conclude that the chances for the bacteria removal in karst channels are very small due to preferential flows and high velocities of the groundwater found inside of them.Mechanisms of bacteria removal from water are conducted by percolation through the surface layer, epikarst and unsaturated zone of karst rock where soil pores and percolation velocities are significantly lower in relation to karst channels.Furthermore, it is determined that intensive precipitation can trigger a major flow of water through the epikarst by virtue of binding smaller fractures and creating a continuous flow.In such conditions, hydraulic pressure can mobilise water stored inside an epikarst towards larger fractures, and eventually saturated zone.Such a flow concept is called piston-type flow [34].Schwarz et al. [34] noticed that the largest number of bacteria reaches the saturated zone shortly after intensive precipitation.The phenomenon is also explained because of smaller fractures binding inside the epikarst and moving of the bacteria sedimented within it.Williams [35] also recognised the piston-type flows when he noticed that the concentration of pollution passing through the epikarst zone exponentially decreases over time.Concentration of the bacteria can be increased with the emergence of "new" precipitation (so-called "pulse" effect).That is why each precipitation can be looked at as a separate event, i.e., as an impulse that triggers pollution towards the saturated zone.

Study Site
The Dinaric karst in Croatia, globally known as "locus typicus" of karst landscapes, is characterised by very irregular karstification which is caused by tectonics, compression, reverse faults and over thrusting structures.The research area is generally located between the Adriatic microplate to the southwest and the Dinaric regional structural unit to the northeast.Their contact is represented by reverse faults that strike northeast-southwest [36].Eocene and Cretaceous limestone is a main component of the terrain [37,38].Carbonate rocks, clastic rocks and quaternary soil are basic rock masses that build the terrain.The most represented rock from those three is the most permeable one, the limestone.The whole hydrogeological system of Bokanjac-Poličnik covers an area of 235.7 km 2 .The Bokanjac-Poličnik aquifer takes around 86% of the Bokanjac-Poličnik catchment, and it spreads around the carbonate rocks formation.The remaining part, around 33 km 2 belongs mainly to non-permeable flysch surfaces [39].
Stratigraphically, the upper zone of the terrain is shallow and highly karstified, while the deeper zone is characterised as a low permeability rock mass, with a network of karst fractures [40].Typical karst shapes in the analysed area are not particularly developed and the most important objects in the hydrogeological sense are active sinks.The most important sink is in Biljane Donje where tracer tests were carried out and where underground connection with the spring Golubinka was proven.It was there that the highest groundwater velocity was measured with apparent values of 8.1 cm/s [41].The flow direction goes from the mid and eastern parts of the catchment towards the shore.The flow begins on the higher parts of the relief, parts of Biljane Donje from where water flows towards northwest, and in that way, it forms the Oko spring.In cases of high waters, the water rises from the Oko spring and forms surface flow called Miljašić Jaruga.The flow then continues towards the Boljkovac and Golubinka springs, and due to "hanging" and complete barriers, part of the water is filtrated south towards Bokanjačko blato and Jezerce well.The tracer tests in the area conducted in 2007 [39] suggests that apparent flow velocities, in the periods of hydrologically low groundwater levels (dry season), are very low compared to other karst terrains.In the periods of low groundwater level, in that morphologically low terrain, gradients are extremely small and water flows extremely slowly.That is why in such dry periods, the groundwater level in wells is rather low, but the inflows are still continuous, and with specific capacity decrease, it can always pump ten to twenty litres per second of the groundwater.Today, for the purpose of water supply, six wells are used.According to the information of municipal corporation in Zadar not a single settlement in the hinterland of Zadar has public sewage system and waste water from the households is treated in septic tanks.Such a way of dealing with waste water presents danger from polluting hydrogeological system Bokanjac-Poličnik and the springs within it.
The data on groundwater levels and concentrations of ENT on the wells in Jezerce and B4 were gathered for the period of 1 January 2012-19 October 2016.Based on [42], an insight was given into the values of the hydraulic conductivity K at the depths from 20 to 50 meters on several locations on the catchment (Figure 1, purple dots).The groundwater level data for that period (1 November 1966-31 June 1967) were also gathered (Figure 1, red dots).

Materials and Methods
The suggested methodology for determining the die-off coefficient of the enterococci in the conditions of transport through the karst aquifer consists of two basic steps (Figure 2).Preliminary analysis, as the initial step of the methodology, serves for making specific assumptions regarding detection of what triggers enlarged ENT concentrations on the analysed wells, as well as regarding potential source of pollution.Information from the preliminary analysis were a basis for making an analytical transport model which is the second step and the backbone of the suggested methodology.The goal of the analytical model was to determine the dominant factor in the process of natural removal of ENT by transportation through a karst catchment.The model consists of two parts of the pollution transport process.The result of the first part (vertical component of the flow) was the value of ENT concentration at the end of the unsaturated zone, C VAD .In the second part, based on the known ENT concentration in the wells, by inverse procedure, the value of ENT concentration at the beginning of horizontal flow part (C SAT ) was calculated.

Time Series Analysis-Transfer Function
Preliminary analysis primarily dealt with the statistical analysis of the gathered data regarding the precipitation, groundwater levels in Jezerce and B4 wells, pumping rates and ENT concentration for the period from 1 January 2012-19 October 2016 at Jezerce and B4 wells.The goal was to determine the correlation strength between the above-mentioned parameters.The model that best describes the relation between input (precipitation, pumping and levels) and output (concentration of enterococci) is defined by the time series analysis.Interpretation of results gave an insight into the dominant variable (precipitation, pumping or levels) which affects the enlarged concentration of ENT on the pumping sites, as well as the approximate reply time of the concentration increase depending on input variables.
A transfer function was used as a method for data analysis, as well as the CAPTAIN Toolbox [43][44][45] integrated in the MATLAB.The transfer function can serve as an effective tool when defining structure and nature of the mathematical models of hydrological systems by using sophisticated statistical methods and interpreting model equations in the field of physics.In statistics, transfer function modelling belongs to the field of time series analysis [44,46] and includes identification and evaluation of model via time series data.The identification of time series data relates to defining the most suitable model structure, i.e., triad values n, m and δ.Parameters n, m and δ can be seen in the Equation (1).Refined instrumental variable algorithms offer a possibility of simple analysis, so-called SISO analyses (one input, one output), but also the complex analysis in a sense of two and more inputs and one output, so-called MISO analyses.Based on the transfer function, it is possible to define time detachment of increased concentrations depending on the initial impulse with a goal to define average time of seepage from potential sources of pollution to the wells.
A simple SISO model can be described by Equation (1) or by a transfer function, Equation (2) [46]: where y t is the output from the linear model at time t, x i is input into the model in time t, δ is the time delay of the input according to the time t, N t is the value of a correlated random variable at time t, e t value of a random variable with standard distribution with zero average and constant scatter at time t.
The user defines the triad n, m and δ, and the algorithm chooses and presents the best possible solutions depending on the chosen identification criterion (YIC, RT2, AIC or EVN).The software offers the possibility of graphic interpretation of results through a graphical user interface where the 20 best solutions are generated according to chosen criteria, and all with the purpose of simplifying identification procedure.
The MISO transfer function, suitable for determining the input parameter with dominant influence on the output, in discrete form is expressed as: where m is the number of input variables that is assumed to influence the output, and A and B i represent polynomials.SISO and MISO analyses were made on the data sets from the Jezerce and B4 wells.

Quantification of the Groundwater Flow Distribution
The numerical model KarstMod was used to estimate the distribution of total seepage through the karst aquifer on the fracture flow (karst channels and fractures) and diffuse flow (matrix).The catchment area of the Golubinka spring was taken into consideration in this step of the preliminary analysis.The quantification of the fracture and diffuse flow component was performed based on the well-known discharge hydrograph from the Golubinka spring and the precipitation data (mm/day) for the period 1 January 2012-31 December 2015.
KarstMod is a platform for the simulation of the precipitation-discharge ratio of the karst springs and the analysis of the hydrodynamic characteristics of the karst aquifer.By defining model structure and its fluxes (precipitation, evapotranspiration, pumping rates in the components), warm-up, calibration and validation periods, an optimal model parameter set can be obtained as the model output.The optimisation parameters of the model and the range of their values are shown in Table 3, together with the adopted optimal values used in the calculation of discharge hydrograph.The two-level structure, with three compartments constitutes the model (Figure 3).Compartment E (epikarst) represents the higher level, while compartments M (matrix) and F (fractures) represent the lower level.It is assumed that the components M and F are filled from the higher-level E and that the main quantities of flow will occur through the fractures (Q FS ).The model has three water balance equations [47]: where E, M and F are water levels in compartments E, M and F respectively, S is the outlet of the system, P e f f is the effective precipitation, assumed as 0.6 • P [48], and Q EM , Q EF , Q MS , Q FS , Q MF are internal discharge rates per unit surface area.Information regarding flow distribution inside the karst aquifer was used to create numerical flow and transport model within software MODFLOW.The numerical strategy SCPE was used for flow simulation.The examination was done here on how numerical model reacts on the distribution of flow onto diffuse flow and fracture flow.Cells which are given relatively high value of hydraulic conductivity represent fractures, while the surrounding cells represent the matrix, and they are given significantly lower values of hydraulic conductivity coefficient.The concept can be visualized as an extremely permeable karst channel (fractures) immersed into a low permeable net called the matrix.For the purpose of checking the flow distribution given by KarstMod, a numerical model of Golubinka catchment was created (Figure 4) with the cell dimensions of 200 × 200 m horizontally, and with two layers vertically.The bottom layer was 1 m thick (−2 to −1 m a.s.l.), and the top layer was 61 m (−1 to 60 m a.s.l.).Slow, diffuse flow took place in the top layer while more significant seepage took place in the bottom layer, where fractures and karst channels were simulated by setting higher values of hydraulic conductivity.Model parametrization in a sense of hydraulic conductivity K x , K y , K z , effective porosity n e f f and total porosity n tot (Table 1) was adopted based on the calibration procedure in [49].
Table 1.Adopted model parameters based on the calibration procedure.

Parameter Upper Layer
Lower Layer

Analytical Model
The analytical model integrated two parts of a pollution transport process (Figure 5); the first part with the percolation from the septic tanks to a saturated zone, and the second part in which horizontal flow and transport were followed inside of the saturated zone.Input data in the implemented analytical transport model were measured groundwater levels on the wells B4 and Jezerce, measured ENT concentration, surface elevation, and data regarding the distance between suspicious pollution sources from the wells Jezerce and B4.For the initial ENT concentration around septic tanks, the value of 4 • 10 8 cfu/100 mL was adopted (this value has been verified by an authorized laboratory).The calculation was conducted through following steps: The calculation was done for each of the suspicious pollution sources and assumed flow directions (Figure 6) with the variation of the parameter-coefficient values in a transport model (die-off and sorption coefficients, turbulent filtration coefficeint, etc.).The transport model parameters/coefficients were varied with the goal to achieve the minimum difference in concentrations ∆C = C VAD − C SAT .Established statistical criterion was the root mean square error (RMSE).MATLAB was used for the implementation of the calculation.

Vertical Flow
The ENT concentration at the end of the vertical flow section C VAD was calculated according to Equation ( 5): where C 0 is the initial ENT concentration (4 • 10 8 cfu/100 mL), T is the on/off function of initial humidity, t p is the percolation time, k tot is the total die-off coefficient, k i is the natural die-off coefficient, k is the sorption coefficient, d is the depth of the unsaturated zone, v VAD is the percolation velocity, and R is the retardation.The on/off function of the initial humidity T considers the impact of the previous precipitation on the wettability of the porous medium, thus ensuring a more pronounced seepage velocity.The parameter T was introduced to cancel the impact of precipitation on the pollution transport after a longer drying season, when subsurface conditions were not suitable for the infiltration and percolation process.The value of the on/off function of the initial humidity was determined according to the following criterion: where P i is the precipitation height in the analysed week, P i−1 is the precipitation height of i − 1 of the analysed week, P i−2 is the precipitation height of i − 2 of the analysed week, P mean is the average precipitation amount for the analysed period of 1 January 2012-19 October 2016.The seepage velocity in the horizontal flow section (v SAT ) was obtained by applying the nonlinear seepage law [50]: where v SAT is seepage velocity, K is the turbulent filtration coefficient, L is the length of the horizontal section of the flow.
The impact of previous precipitation on the value of the ENT concentration is considered through modification of the Equation ( 5), so the final ENT concentration at the end of an unsaturated zone C VAD,i was calculated by using the following equation: where C VAD,i is the ENT concentration at the end of the unsaturated zone, C VAD,(i−1) is the ENT concentration at the end of the unsaturated zone for i − 1 of the analysed week, C VAD,(i−2) is the ENT concentration at the end of the unsaturated zone for i − 2 of the analysed week, F 1 and F 2 are the coefficients of the delayed concentration.
The second and the third member on the right side of Equation ( 10) refer to the delayed concentration, i.e., the one that is left behind in unsaturated zone.They depend on catchment conditions, i.e., seepage velocity and travel time from the pollution source to the wells B4 or Jezerce.F 1 and F 2 are also calibrated and optimised parameters with the possible values expressed hereafter: where t S is the seepage time in the saturated zone.

Horizontal Flow
The aim of this part of the model is to calculate ENT concentration at the beginning of the horizontal part of the flow C SAT (Equation ( 13)), based on the measured values of concentration on the wells Jezerce and B4.
where C B4 or C JEZ are measured ENT concentrations on wells B4 and Jezerce.The value of the k i is equal to the value of the die-off coefficient in the Equation ( 6).Seepage time in the saturated section of the flow is calculated using the following equation: It should be mentioned that the calculation was applied only in cases when the condition of total time of seepage through the karst medium was satisfied in the following way: t s + t p < 7 days (time from the initial rain impulse to the end of flow through horizontal part of the flow is less than seven days).If it is the opposite way, the concentration in an observed time is zero and it is eliminated from further calculation.
The time limitation was established because of maximal correlation strength accomplished while adopting seven days lag between precipitation occurrence and increased ENT concentration on wells (from the time series analysis).Delayed concentration in the unsaturated zone were considered through adequate members in the Equation (10).

Time Series Analysis-Transfer Function
Three SISO analyses (precipitation/groundwater levels/pumping rate-ENT concentrations) and one MISO analysis (precipitation/ groundwater levels-ENT concentrations) for Jezerce and B4 were made.From the results of the SISO analysis (Figure 7) it is evident that precipitation have the greatest impact on the increased ENT concentration in the wells.Optimal values of n, m and δ as well as the value of the identification factor are shown in Table 2.The results represent the best compromise between the identification criterion RT2 and YIC.The analysis shows that the most pronounced correlation strength between precipitation and concentration is obtained for the time lag δ = 5 days for the B4 and seven days for the Jezerce well.This value corresponded to the velocity of the groundwater flow of 1-3 cm/s from the Poljica, Poličnik, Visočane, Briševo settlements to the wells B4 and Jezerce, obtained with the tracer tests [40].These settlements were taken into consideration in the further implementation of the methodology.In Figure 8, the results set of SISO and MISO analysis is shown, together with the measured ENT concentrations.The calculated discharge hydrograph and the results of flow distribution calculation, obtained based on the optimal parameters from Table 3, are shown in Figures 9 and 10.According to the results of the model, 19% of the total flow took place in the diffuse phase (matrix) and 81% in the fracture component of the flow (fractures and karst channels).
For the control of the flow distribution in MODFLOW, the cross section marked in Figure 4 was taken into consideration.Yellow cells represent the matrix (K x , K y , K z = 1 • 10 −4 ), and blue cells represent karst channels, i.e., fractures (K x , K y , K z = 3 • 10 −1 ).In cross-sectional view, there were 14 cells altogether, 10 of which simulated the matrix.The average velocity inside of the matrix was 5 • 10 −5 m/s.Total flow per unit meter was estimated through the product of velocity and height of the water level inside of the cells and amounted 0.03 m 3 /s for seeping through the matrix.In the fractures (blue cells in Figure 4) the average velocity amounted to 0.042 m/s.Total flow per unit volume through fractures was 0.17 m 3 /s (velocity times cell height of 1 meter times number of cells).The given relation of total flow apportionment was 84% for seeping through channels and 16% of water which was seeped through low permeable matrix.Those results are in accordance with the results from KarstMod, and therefore the SCPE strategy can be applied in creating the numerical flow and transport model for the Bokanjac-Poličnik catchment.

Analytical Model
The value range of variable and fixed parameters that appear in Equations ( 5)-( 14) were firstly defined (Table 4).Flow directions from Figure 4 (distances) and the initial ENT concentration around septic tanks (4 • 10 8 cfu/100 mL) were fixed values.ENT concentrations in the Jezerce and B4 wells were obtained by measurement.Die-off coefficient k i depends on water temperature and varies from 0.15 day −1 (10 °C) to 0.5 −1 (20 °C) [27].That was consistent with water temperature of the analysed aquifer over the year so the die-off coefficient k i in the calculation was taken within the limits of those values, which means 0.15 in the winter period, 0.5 in the summer period, and linearly interpolated in between.Retardation coefficient for faecal indicator bacteria acquired values in a range from 3 to 10 [51], and within the analytical model represented the optimised parameter.When it comes to percolation velocity v VAD , [35] recognised three phases of flow in the unsaturated zone; zone of rapid velocities from 0.5-2 cm/s, second phase with velocities from 10 −2 cm/s and the third with velocities lower than 10 −3 cm/s.Percolation velocity v VAD was also an optimised parameter in the calculation.Sensitivity analysis on the statistical criterion RMSE for the variable ∆C was suggested because literature data of some parameters had a relatively wide value range.Such implementation was suggested with the aim to narrow the value range of certain parameters, and preferably, to fix the least sensible.After that, Monte-Carlo simulation was used to generate 3 • 10 5 combination of values of variable parameters, and then they were implemented in Equations ( 5)- (14).Then, the calculation ∆C = C VAD − C SAT was made within the simulation period of 1 January 2012-19 October 2016.Using the RMSE statistical criterion for variable ∆C, the optimal combination of the parameter values for each path 1-8 (Figure 6) was obtained (Tables 6 and 7). Figure 13 shows the calculated concentration of C VAD and C SAT for the paths 1 and 4 in the direction of the Jezerce well, gained with the adopted optimal parameters from Table 6.The Figure 13 shows that C VAD and C SAT concentrations are of the same magnitude and what is important, the pulse of the concentration occurrence appears to be the same for both phases of the transport, i.e., there was no delay in the concentration appearance in relation to each other.6).

Discussion and Conclusions
This paper presents a methodology for determining the die-off coefficient of faecal indicator bacteria, ENT, in terms of transport through the karst environment.The methodology assumes the successive application of statistical, numerical and analytical modelling in order to define the self-purification capacity of the analysed aquifer (the effects of sorption and die-off expressed by the corresponding coefficients).
With the transfer function, used as a method of time series analysis, a specific answer to the question of a dominant trigger of the enlarged ENT concentrations in the wells was obtained.In addition, time series analysis gave an insight into the average groundwater velocities later used to determine the potential source of pollution.In the applied analytical model, the calculation of ENT concentration removal was divided into two parts, the unsaturated and saturated transport phase.Removal of ENT concentrations in the unsaturated zone, due to the low percolation velocity and small size of the fractures, binds to two processes: natural die-off and sorption.By transporting through a saturated environment where much higher seepage velocities occur, the impact of sorption is negligible and the dominant process of ENT removal is a natural die-off.
The influence of the individual parameters used in the analytical model was analysed by sensitivity analysis.The diagram in Figure 11 shows the influence of the turbulent filtration coefficient value at variable ∆C.It is evident that ∆C is smaller in the case of a turbulent filtration coefficient of K = 0.5 m/s, which is also confirmed by the lower values of the statistical criterion RMSE.Also, the calculated average seepage velocity (Equation ( 9)) for the conductivity coefficient K = 0.5 m/s for the time period 1 January 2012-19 October 2016 is 1.26 cm/s and maximum 2.3 cm/s.Such values are consistent with the apparent velocity values in the area obtained from tracer tests [40].Therefore, in further analysis, the value of the turbulent filtration coefficient K = 0.5 m/s was used.The sensitivity analysis for the variable v VAD was conducted for the range of 10 −4 to 10 −2 m/s.From Figure 12 the exponential increase of the RMSE is clearly visible from the value v VAD = 10 −3 m/s upwards.The difference of the statistical criterion RMSE between 10 −3 m/s and 10 −2 m/s is approximately seven orders of magnitude (RMSE ≈ 6.5 • 10 7 for v VAD = 10 −2 m/s, RMSE ≈ 20 for v VAD = 10 −3 m/s).The analysis showed that the percolation velocity of 10 −4 m/s results in the complete removal of faecal indicator bacteria to the contact with the saturated zone due to the extremely slow percolation which leaves enough time for ENT removal.With such velocity values, ENT to the saturated zone can only break through if ignoring the sorption effect, i.e., if k = 0 (Equation ( 6)).On the other hand, due to the velocity of 10 −2 m/s, the percolation is rapid and ENT encounter the saturated zone in virtually unchanged concentration.Therefore, value of the percolation velocity of v VAD = 9 • 10 −4 m/s was adopted, which, according to [35], corresponds to the mean velocity in three recognised phases of percolation.
The mean value of the total die-off coefficient by transport through the unsaturated zone in the analysed case is k tot = 8.25.Within the saturated zone the total die-off coefficient k tot is within the limits of 0.1 and 0.5.From the results of the analysis it can be concluded that the unsaturated zone represents a suitable environment for ENT removal (reduction of six orders of magnitude in relation to the initial concentration) where sorption is the dominant process of natural ENT removal.
Further research is planned to focus on the unsaturated zone and the processes within it, especially from the aspect of removing pollution from water.The significance of the unsaturated zone in pollution removal from the water has been proved by the implementation of this methodology.

Figure 1 .
Figure 1.Schematic hydrogeological map of the Bokanjac-Poličnik catchment area with the locations of boreholes, piezometers, wells and groundwater barriers.

Figure 2 .
Figure 2. Basic steps of the proposed methodology for determining the die-off coefficient (k tot -total die-off coefficient, k i -natural die-off coefficient, k-sorption coefficient).

Figure 3 .
Figure 3. Structure of the KarstMod model platform.The 90-day period (1 January 2012-31 March 2012) was used as the warm-up period and these results were not used for the model calibration.For the calibration, 1004 days were used (1 April 2012-31 December 2014) for which the optimum model parameter set was obtained.For the model validation, the year 2015 was used.Model performance was tested by using the Nash-Sutcliffe efficiency coefficient (NSE) and the modified Balance Error (BE).Information regarding flow distribution inside the karst aquifer was used to create numerical flow and transport model within software MODFLOW.The numerical strategy SCPE was used for flow simulation.The examination was done here on how numerical model reacts on the distribution of flow onto diffuse flow and fracture flow.Cells which are given relatively high value of hydraulic conductivity represent fractures, while the surrounding cells represent the matrix, and they are given significantly lower values of hydraulic conductivity coefficient.The concept can be visualized as an extremely permeable karst channel (fractures) immersed into a low permeable net called the matrix.For the purpose of checking the flow distribution given by KarstMod, a numerical model of Golubinka catchment was created (Figure4) with the cell dimensions of 200 × 200 m horizontally, and with two layers vertically.The bottom layer was 1 m thick (−2 to −1 m a.s.l.), and the top layer was 61 m (−1 to 60 m a.s.l.).Slow, diffuse flow took place in the top layer while more significant seepage took place in the bottom layer, where fractures and karst channels were simulated by setting higher values of hydraulic conductivity.Model parametrization in a sense of hydraulic conductivity K x , K y , K z , effective porosity n e f f and total porosity n tot (Table1) was adopted based on the calibration procedure in[49].

Figure 4 .
Figure 4. Spatial distribution of hydraulic conductivity coefficients K x , K y , K z in the bottom layer (from −2 to −1 m a.s.l.) with the marked section for flow distribution control and boundary conditions (Golubinka and P1).Blue area indicates the values K x , K y , K z = 3.0 • 10 −1 , and the white area indicates the values K x , K y , K z = 1.0 • 10 −4 .

Figure 5 .
Figure 5.The scheme of the analytical transport model.
unsaturated zone (vertical flow)-calculation of the ENT concentration at the end of the unsaturated zone C VAD (the place of shift from vertical to horizontal flow); (b) Transport through the saturated zone (horizontal flow)-calculation of the ENT concentration C SAT at the beginning of the horizontal flow route (inverse procedure) based on measured concentration C B4 and C JEZ on the wells B4 and Jezerce.

Figure 6 .
Figure 6.The location of the pollution source in relation to the wells Jezerce and B4 along with possible flow paths.

Figure 7 .
Figure 7.One input, one output simple analysis (SISO) analysis of the input data and enterococci (ENT) for wells B4 and Jezerce.

Figure 8 .
Figure 8. Results of SISO and two and more inputs, one output complex analysis (MISO) for the Jezerce well.

Figure 9 .
Figure 9. Discharge hydrograph of the Golubinka spring obtained by the KarstMod model simulation.

Figure 10 .
Figure 10.Flow distribution between the diffuse (MS) and the fracture (FS) flow components through the analysed karst aquifer.

Figure 13 .
Figure 13.Diagrams of C VAD and C SAT concentrations for flow paths 1 and 4 towards Jezerce well obtained with optimal combinations of the variation parameters values (Table6).

Table 2 .
Model structure and values of the identification factors.

Table 3 .
Ranges of parameter values carried out by the optimisation procedure based on the discharge hydrograph of Golubinka spring and the optimal values of the parameters used for calculating hydrograph with the KarstMod model

Table 5 .
Values of variational and fixed parameters after the sensitivity analysis.