Abstract
This paper shows the results of the hydraulic-hydrologic calculations of karst spring discharges and the groundwater level in the aquifer of spring catchment. The calculations were performed for the Golubinka spring in Zadar area for the 4-year period. The chosen approach was a model using relatively scarce data set, including limnigraphic data on the difference between the sea water level and the freshwater level on the spring itself and the precipitation data from the meteorological station Zadar. The determination of discharge hydrographs, based on inherent assumptions and available data, yields the proportionality coefficients between the discharge and the limnigraphic data on the Golubinka spring itself. Further, based on the discharge hydrograph, groundwater level oscillation was determined. The resulting spring discharge hydrograph and groundwater levels, along with the assumption of Golubinka spring as the only spring on the catchment, were used in creating turbulent seepage model of the fractured system within the aquifer, which evidently extends along the axis of the Golubinka spring catchment. The model yielded suitable turbulent seepage coefficients of the fracture system. By using the numerical model KarstMod it was estimated that, on average, concentrated fracture flow drains around 85% of infiltrated volumes and the remaining 15% accounts for diffuse matrix flow. Finally, the Modflow model was used in order to get insight into the flow field and the distribution of equipotentials in the aquifer of the Golubinka catchment.
1. Introduction
The Dinaric karst in Croatia, globally known as “locus typicus” of karst landscapes, is characterized by very irregular karstification which is caused by tectonics, compression, reverse faults, and overthrusting structures. Due to the morphological and hydrogeological differences, Dinaric karst itself can be distinguished into many types. Segregation is necessary due to heterogeneity of the rock mass and nonuniform infiltration that varies over space and time. Having that in mind favors a more detailed and precise analysis of karst terrains, when dealing with flow rates estimations, contaminant monitoring and tracking, or creating a protection model for karst aquifer [1]. Because karst in the Ravni Kotari area near Zadar city differs from karst in the Lika region and karst on the Adriatic islands, flow model needs to be adjusted to the natural conditions [2] and should not be generalized over whole area of Dinaric karst. In order to identify recharge areas and infiltration processes in karst, the geomorphological and the topographical surveys are needed. Independently of soil covers, areas of karstified carbonate rocks are classified as zones of autogenic recharge. For example, 80% of precipitation can infiltrate bare karstified limestones, that is, karrenfelds [3]. The increase of vegetation cover significantly increases evapotranspiration from karst terrain so the amount of recharge decreases. In addition, climatic conditions such as the intensity, type, and amount of precipitation and temperature also affect the percentage of recharge. Knowing the hydrodispersion characteristics of the aquifer, boundary conditions, hydraulic heads, and spring and sinks locations is necessary to acquire insight into existing conditions on the catchment. It is expensive and complicated to collect all the necessary data, which represent a problem when dealing with a numerical analysis, where these kinds of data are very needed. However, the problem can be reduced by using appropriate assumptions and relying on the researcher’s experience and previous knowledge.
The scale of observations can significantly influence the opinion about the aquifer and its characteristics. For example, hydraulic conductivity value is especially sensitive to the scale effect [4], which should be taken into account due to importance of this coefficient in assessing the rate of groundwater flow. The aquifer hydraulic conductivity determined by the pumping tests can indicate more than five orders of magnitude greater value than the hydraulic conductivity value of the isolated matrix structure [5]. Matrix structure has a small primary permeability, due to pores which are usually smaller than 10 micrometers. Even highly porous soil, for example, chalk and some limestone, is not very permeable, and its specific retention is high. Because of such small permeability, the flow through a matrix structure is characterized as a diffuse flow, accounting for a minor part of the total groundwater flow. The largest part of the groundwater flow in the karst environment is conducted through fractures, where higher rates of groundwater velocities are developed. As reported by Milanovic [6], in the classic Dinaric karst of Herzegovina, the most frequent velocity is between 864 and 1728 m/d, based on the 281 tests of dye tracing. These values represent apparent velocities developed in the fractures of the aquifer. In order to assess the groundwater flow directions and the overall flow characteristics in a real karst aquifer it would be insufficient to use only the hydraulic head information. The measurements of hydraulic head should be combined with field tests, such as dye tracing, and the overall understanding of various hydraulic factors, such as flow through pipes and “nonhydraulic” factors such as geomorphology, limestone (carbonate) sedimentology, tectonics, and geologic history, including paleokarstification. Darcy’s equation, which is applicable for intergranular porosity and laminar conditions, is not suitable for fracture porosity. In karst, flow takes place through micro and macro fractures, through interconnected pores within a matrix structure, but also through solutional cavities of any imaginable shape [3]. Therefore, it is quite difficult to determine the “representative” parameter values for the Darcy equation; the hydraulic conductivity , the hydraulic gradient , and the entire cross-sectional area of flow , inclusive of voids and solids. In order to get more realistic results of the groundwater flow rates, it is necessary to integrate the equations of turbulent flow through various types of fractures, pipes, and channels and Darcy’s equation for laminar flow through the rock matrix.
The aim of this paper is to model the flow through the Golubinka spring catchment despite the scarcity of measured data. The idea is to get acquainted with aquifer characteristic and conditions within (flow rates and distribution, velocities, hydraulic conductivity coefficients, etc.), so this knowledge can be applied in further researches.
2. Methodology
In this research we were forced to deal with a lack of information about the analyzed catchment area. The placement of the fractures and preferential flow paths were not fully obtained or known. There was no information available considering the underground water levels, with exception of the field research provided by the IGH institute in 2013, when exploratory drilling was performed on the analyzed catchment on 29.9.2013. The main goals of this research are to obtain a discharge hydrograph for karst spring Golubinka, to determine a groundwater level hydrograph and seepage quantities within aquifer of Golubinka catchment area, using inherent assumptions and hydraulic equations. The research is designed to be successive: the results of mentioned calculations will serve as an input for the two numerical models, one regarding quantification of flow components and the other used for simulation of flow fields. Available data used as the main input for this research were the precipitation data obtained from a nearby meteorological station and data from the limnigraphic station located at the Golubinka spring. The analysis was performed for period between 1.1.2012 and 31.12.2015. Since there is no data on the Golubinka spring discharge, the spring discharge hydrograph is synthesized with limnigraphic data on the difference between the sea water level and the freshwater level on the spring itself. It was assumed that the total infiltrated volume of precipitation equals the total volume of discharge at the Golubnika spring. This allowed us to analyze this area as a closed system, where Golubnika spring is the only place of discharge. For the calculation of the discharge hydrograph, two hypotheses were tested, one which assumes a linear correlation between the discharge and the water level difference, and the other which assumes a nonlinear correlation.
The oscillation of groundwater level in aquifer of the Golubinka catchment is simulated, on the basis of the continuity equation proposed by Katsanou et al. [7], assuming infiltration and storage coefficients and using the spring discharge values obtained with the linear and nonlinear approaches. There was only one in situ data (piezometer ; IGH, 2013) considering the water level at the analyzed period, used as a control point in the calculations. The series of groundwater level measurements from 1967, which we had at our disposal, were used for comparison of the obtained results. Furthermore, turbulent component of overall seepage through the karst aquifer of the Golubinka spring catchment area has been simulated, using the results from previous calculations, that is, spring discharge quantities and groundwater levels. Nonlinear Manning’s equation was used for the calculation of turbulent seepage through the fracture system of the Golubinka spring catchment. Seepage hydrograph of the aquifer will be compared with the Golubinka spring discharge hydrograph obtained with the first step of the methodology, in order to find the suitable turbulent seepage coefficient for this fracture system. In that way, Golubinka spring hydrograph will serve for calibration of seepage hydrograph and turbulent seepage coefficient () will be the main outcome of this calculation. This coefficient will be used later in numerical modeling of the analyzed aquifer. It is expected that the majority of flow will occur in the fracture section of the aquifer [8–10], but the question is in what proportion will this happen compared to the diffuse matrix flow. The numerical model KarstMod is utilized for estimation of distribution between the fracture component and the matrix component of overall groundwater flow, using the precipitation time series, pumping rates from the Golubinka spring, and the discharge hydrograph obtained by the linear and nonlinear approach. Finally, the Modflow model will be used in order to get insight into the flow field and the distribution of equipotentials in the aquifer of the Golubinka catchment.
3. Study Site
The Golubinka spring, located in the northern Dalmatia in Croatia, is a typical karst spring (Figure 1). The catchment is mostly composed of limestones from the Cretaceous and Palaeogene era with some dolomitic parts. The hydrogeological barriers are usually “hanging” barriers composed of Quaternary sediments or flysch and below which the groundwater flow is formed. Near the Golubinka spring there is a typical fracture flow which is well connected with the periphery of the catchment, which was proved with tracer tests. 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 [11, 12]. Eocene and Cretaceous limestone are a main component of the terrain [13, 14]. Miljašič Jaruga is the only temporary surface stream in the Zadar hinterland. The catchment area of the Golubinka spring (Figure 1) is an integral part of the wider catchment area that covers the entire hinterland of the city of Zadar, called Bokanjac-Poličnik. The most known karst springs and pits, some of which are used for water supply, are located in this area. Hydrogeologically, the catchment of the Golubinka spring can be separated from the rest of the Bokanjac-Poličnik catchment by zonal underground water divide, which is unfortunately known only approximately (Figure 1). In the coastline area, the carbonate rocks are in contact with sea water in quite narrow area, because flysch rock in Ljubački bay acts as an undersea barrier [15]. In that way the groundwater flow can suppress intruding sea water until the end of the summer period when groundwater level and discharge drop. In the areas formed by carbonate rock the effective infiltration is very large and in the areas of impervious terrain the effective precipitation ends in the aquifer through sinks [16]. The intermittent surface flows that arise on the impermeable flysch rock also end up in sinks. The Biljane Donje sink is located in the periphery of the Golubinka catchment (Figure 1). Another groundwater barrier is located southwest of the Golubinka spring, where there is a five-meter thick layer of residual soils [17]. However, some water flow passes under residual soil layer towards the Bokanjac spring. Stratigraphically, the upper zone of the terrain is shallow and highly karstified, while the deeper zone is characterized as a low permeability rock mass, with a network of karst fractures [16]. The observed direction of the groundwater flow in most cases is northwest. This is an underground connection between the Biljanje Donje sink and the Golubinka spring with the groundwater velocity of 8.1 cm/s, which has been demonstrated through tracer tests [16]. The measured groundwater level in piezometers by the Oko spring (Figure 1) is around 60 m a.s.l. That causes a significant hydraulic gradient and high groundwater velocities towards the northern part of the catchment area, that is, the Golubinka spring or Boljkovac spring. The measurements on the study area from 1961 to 1990 [18] indicate a positive correlation between the average monthly precipitation and the discharge at the Golubinka spring. The highest monthly discharges occurred in February, with an average of 1.000 l/s, and the lowest were in July, with an average of 100 l/s [18]. The average annual discharge of the Golubinka spring was 417 l/s and the average annual precipitation was 992 mm/year.
4. Determining the Discharge Hydrograph for the Golubinka Spring and the Groundwater Level in Aquifer
In the analyzed period (1.1.2012–31.12.2015) the daily limnigraphic data and the daily precipitation data were available (Figure 2). The limnigraphic data are the water level differences between the “sea” side and “freshwater” side of the partition wall built on the spring itself. The precipitation data are from the meteorological station Zadar (Figure 1). It should be noted that these data do not relate to the overflow height over the calibrated flume, but only to the fresh water surface elevation above the sea level on the place of the spring next to the sea. Consequently, there is no reliable hydrograph for the Golubinka spring.
To obtain the discharge hydrograph we introduced the coefficient of the Golubinka spring discharge and the water level difference from the “sea” and “freshwater” side of the partition wall . The coefficient depends on the estimated nature of their relationship. The two hypotheses were tested. The first one assumes linear correlation between spring discharge and water level difference and the second assumes nonlinear correlation . The values were obtained based on the assumption that a linear correlation between the precipitation and the average annual discharge from the spring exists. Data from 1961–1990 [18] assist in obtaining spring discharge values. From that period we adopted the average annual precipitation of 992 mm/year and the average annual spring discharge of 417 l/s and determined an average discharge for the analyzed period (2012–2015) on the basis of known data on average annual rainfall for the analyzed period (1036 mm/year). The simulation was performed until the same precipitation and discharge ratio of the analyzed period and the one from [18] were reached (). The calculated hydrographs for utilizing linear (for ) and nonlinear (for ) correlation are shown in Figure 3. The average discharge for both cases in the analyzed period is 435 l/s, but the peak values differ. The peak discharge at the spring for the period 2012–2015 is 1364 l/s for linear correlation and 831 l/s for the nonlinear correlation.
In order to determine the groundwater levels for the same period a simple approach based on the continuity equation (1) was used [7]. Accordingly, it is assumed that the groundwater level on the reference piezometer position (Figure 1), located in the middle of the catchment, is the reference value for the entire catchment area.where , are groundwater level at the beginning and the end of the calculation [L], is discharge at the Golubinka spring [L3T−1], is pumping rate [L3T−1], is infiltration coefficient (adopted 0.6; [19]), is catchment area (adopted 65 km2; [18]), and is storage coefficient (adopted 0.01; [20]).
The infiltration coefficient is usually in the range 0.3 to 0.9 and the characteristic value for the Dinaric karst is around 0.6 [20]. At our disposal we had measurements of groundwater level from 1967 [21] which assist as a comparison for the results obtained for the analyzed period. The mean groundwater level in piezometer in 1967 was 53 m a.s.l. and the maximum was 64 m a.s.l [21]. The total precipitation in 1967 was 1762 mm which is significantly higher than the annual average precipitation of 992 mm in the period 1961 to 1990 [18] and of 1036 mm in the analyzed period 2012 to 2015. Therefore, the mean and the maximum groundwater levels on the reference piezometer in 1967 should be higher than those calculated for the period 2012 to 2015. In addition, it is necessary to adopt an initial groundwater level for the calculation. It should be noted that during the simulation period 1.1.2012 to 31.12.2015 there were no continuous measurements of water levels in any piezometer within the analyzed catchment, except for piezometer ( = 15019′47,8′′E; = 44011′19,2′′N; altitude of 55 m, Figure 1) where, during the starting period of the hydraulic conductivity test pumping on 29.9.2013, the registered groundwater level was 49.5 m a.s.l. [22]. Therefore, the initial condition of groundwater level is adopted in a way that the calculation of groundwater level gives the measured value of 49.5 m a.s.l. on 29.9.2013 at the location of piezometer . Figure 4 shows the calculated, (1), groundwater levels for the analyzed period 2012 to 2015 with adopted model parameters: = 65 km2, = 0.01, and = 0.3. The calculation of groundwater level obtained on the basis of discharge hydrograph with assumed linear correlation, , yields the peak groundwater level of m a.s.l. and the mean level of m a.s.l. For the nonlinear correlation the calculations yield higher peak level of = 62.35 m a.s.l. and higher mean levels of m a.s.l. As previously stated, calculated values of and for both linear and nonlinear assumption are lower than the ones measured in 1967 (Table 1). The groundwater measurements on piezometer in 1967 during strong precipitation events of 374 mm in the seven-day period (4.11.1967 to 10.11.1967) showed a 12 m increase of groundwater level. Through the seven days of the analyzed period (1.9.2012–7.9.2012) when total precipitation was 106.3 mm, increase of groundwater level in case of the linear model was 3.17 m, and 3.1 m in case of the nonlinear model. These values correspond to the measured increase of groundwater level over seven days in 1967. The structure of (1) indicates the possibility of obtaining the same resultant groundwater level hydrograph (Figure 4) if the infiltration and the storage coefficient are multiplied by the same number the area of the basin is divided by (i.e., for the number 2: ; ; km2).
5. Seepage Model Adoption for the Karstified Aquifer of the Golubinka Spring Catchment Area
The following paragraph describes the underground water seepage conditions through the catchment area of the Golubinka spring. With specific seepage of karstified aquifers, it should be noted that for seepage the linear form of Darcy’s law is applicable for cases with critical Reynolds number below 1. For conditions where Reynolds number is above 1 a transitional and turbulent flow occurs in which resistance is nonlinear so the slope of piezometric height is expressed in the form . For the turbulent flow and the appropriate seepage law is , where is the turbulent seepage coefficient. The groundwater flow in the catchment area of the Golubinka spring occurs through an uneven fracture system with different conductivity values along the flow axis. For the calculation of the turbulent flow in the fractures, the Darcy-Weisbach or the Manning equation can be used: , where is the relative cumulative function of conductivity, in units of m3/s. If it is imagined that there are one or more natural underground fracture systems in cross section of aquifer, then the flow is calculated from the general cumulative conductivity and the slope of the piezometric head. In an uneven fracture system conductivity varies along the axis of the flow, wherein the mean conductivity, (2), equals the arithmetic mean of the conductivity at the edges of the section: and . was calculated using expression: , where is cross section of the aquifer, . By integrating the above-mentioned expression for turbulent flow in the fractures, and by separating the variables, (3) was obtained.where is the arithmetic mean of the conductivity [L3T−1], and are hydraulic conductivities [L3T−1], is distance between the edges of the section [L], is groundwater level at piezometer [L], is groundwater level at Golubinka spring [L], and is turbulent seepage within aquifer [L3T−1].
Since we do not have any discharge data from other springs in Ljubački bay we have adopted an assumption of Golubinka spring as the only place of discharge of the catchment. That, along with assumption of equality of total infiltrated and discharged volume of water, allowed us to calculate aquifer seepage hydrograph on the basis of spring discharge hydrograph. By using (2) and (3) the seepage through the fracture system of the Golubinka catchment area is calculated with the variation of (Figure 5), which represents the calibration parameter. The following values were used: = 10 km (distance to Golubinka spring, Figure 1), = 2000 m (the width of the fracture system at the piezometer ), = 150 m (width of the fracture system at the Golubinka spring), is of Figure 4, and is of Figure 2. The results indicate a significant level of correlation between the groundwater flow and the discharge from the Golubinka spring, Figure 5. For the assumption of a linear model the calculated conductivity was m/s (at piezometer ) and m/s (at Golubinka spring), and for the assumption of a nonlinear model the results were m/s and m/s.
(a)
(b)
The question of the distribution of flow between the diffuse matrix and the fracture component of the flow through the karst aquifer remains. The answer to this was given by analyzing it through the numerical model KarstMod.
6. Cascade Reservoir Model for Quantification of the Groundwater Flow Distribution
The KarstMod model is used for the calculation of rainfall-discharge relationship of the karst springs and for the hydrodynamic analysis of individual compartments and their interaction in wider model. By defining model structure and its fluxes (precipitation, evapotranspiration, and pumping rates in the components), warm-up, and calibration and validation periods, an optimal model parameter set can be obtained as the model output. The two-level structure, with three compartments, constitutes the model (Figure 6). Compartment (Epikarst) represents the higher level, while compartments (Matrix) and (Fractures) represent the lower level. It is assumed that the components and are filled from the higher level and that the main quantities of flow will occur through the fractures (). The model has 3 water balance equations [23]:where , , and are water levels in compartments , , and , respectively, is outlet of the system, is the effective precipitation, assumed as [19], and , , , , and are internal discharge rates per unit surface area.
Discharge at outlet is defined aswhere is the recharge area of the catchment and is the discharge rate abstracted at the outlet per unit surface area.
Exchange function is defined aswhere is the specific discharge coefficient, is a positive exponent, is the algebraic flow from compartment to compartment , and is a unit length.
Input data of the model were precipitation, pumping rates from the Golubinka spring, and the discharge hydrograph from Figure 3. The 90-day period (1.1.2012–31.3.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.4.2012–31.12.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 ().where stands for either the discharge or the square root of the discharge and represents observed data.
Both coefficients, NSE and BE, can range from to 1. The greater values of NSE and BE coefficient correspond to a better match between the model and observations. The multiobjective function, WOBJ, defined as the weighted sum of two of the performance criteria , was used for the performance testing during the model calibration. Monte-Carlo simulation with a Sobol sequence sampling of the parameter space was used in the calibration procedure, which took place until all parameter sets that satisfied the criteria are collected. Since more sets of parameters satisfied the criteria during the calibrations, the additional evaluation was used. The KarstMod model evaluates the uncertainty of the model results by the approach derived from the Regional Sensitivity Analysis (RSA) and the Generalized Likelihood Uncertainty Estimation (GLUE). In this way each set of parameters with are used in the prediction process, where the WOBJ value is used as a likelihood measure for each behavioural parameter set [23]. The KarstMod then proposes simulation results of discharge at time , with 90% confidence interval, as shown in Figure 9. The value range for the parameter optimization is given in Table 2, together with the adopted final “optimum” values used in the calculations and for which the obtained results are shown in Figures 7 and 8.
(a)
(b)
(a)
(b)
7. Creation of the Hydrodynamic Model with the Porous-Equivalent Approach
In order to gain insight into the flow field and the spatial distribution of equipotentials for the Golubinka catchment aquifer, the deterministic model Modflow was used [24], which is often used for the simulation of groundwater flow in various hydrogeological conditions [25, 26]. Various types of deterministic models have been applied to karst aquifers. The simplest approach is often referred to as the porous-equivalent media model approach (also called the single-continuum porous-equivalent or SCPE). Simulations of advective transport using a SCPE model are infrequently performed with varying degrees of success [27–32]. In these studies, groundwater Darcy flux was mainly governed by the model cell hydraulic conductivity and total model layer thickness. The effective porosity is used for calculation of a pore velocity that matches times of travel from tracer tests. In some cases where conduit locations were known, the finite-difference cells with conduits were assigned much greater hydraulic conductivity values than surrounding cells and were successful in reproducing transient spring discharge and matching tracer test travel times.
The potential flow equation, used in the basic version of Modflow (SCPE model), assumes that flow is laminar and that temperature, density, and viscosity are constant over the model domain.where , , and are hydraulic conductivity coefficients in the , , and coordinate axes, is piezometric level, is flow (volumetric) per unit of volume that represents the sinks and/or springs, the is coefficient of the specific storage in the porous media, and is time.
The governing equation for the given initial and boundary conditions is solved with the finite-difference method. The Golubinka catchment was presented by an orthogonal network with grid spacing of 200 × 200 m (Figure 10) in the horizontal plane with three sublayers in vertical direction. The bottom layer has a thickness of 1 m (from −2 to −1 m a.s.l.), the middle layer 51 m (−1 to 50 m a.s.l.), and the surface layer 15 m (50 to 65 m a.s.l.). The slow, diffuse flow takes place in the middle and the surface layer, while the much more pronounced seepage appears in the bottom layer. By using the initial groundwater level at 30.4 m a.s.l. the flows were calculated for the period 1.9.2012–26.05.2015 (1000 days) by using the boundary conditions of overlay thickness variation of fresh water (Figure 2) at the Golubinka spring, for the same daily precipitation amounts and the same infiltration coefficient of 0.6.
Parameterisation of the model in terms of determining appropriate values of hydraulic conductivity coefficients , the effective porosity , and total porosity coefficient is adopted on the basis of the calibration procedures which rely on groundwater level at piezometers (Figure 4). The adopted values are shown in Table 3 and Figure 11. The results for the groundwater level comparison for piezometer , between calculated and modeled values, are shown in Figure 12. Figure 13 shows the flow pattern for the fracture (bottom) layer of a thickness of 1 m (from −2 to −1 m a.s.l.) for 15.11.2014, at the onset of the greatest groundwater level at the position of piezometer , during the simulation period 1.9.2012–05.26.2015.
8. Conclusion
This paper shows the results of the hydraulic-hydrologic calculations of Golubinka karst spring discharges and the groundwater level in the aquifer of spring catchment for the 4-year period. The chosen approach was a model using limnigraphic data on the difference between the sea water level and the freshwater level on the spring itself and the rainfall data from the meteorological station Zadar. Two hypotheses were tested, one which assumes a linear correlation between the discharge and the water level difference, and the other that assumes a nonlinear correlation.
Calculation procedures yield the proportionality coefficients between the flow and the limnigraphic data on the Golubinka spring itself. Obtained discharge hydrograph along with the groundwater level of referenced piezometer assists in construction of seepage model for karstified aquifer of the Golubinka spring catchment area. In the model, a turbulent flow was assumed and Manning equation for flow was used. As a result of the simulation the turbulent seepage coefficient was gained. The linear model shows a more significant level of correlation with seepage model, compared to nonlinear one. It can be concluded that the linear relationship between discharge and height is more suitable for determining the spring discharge and, later, seepage through the aquifer. The average error and root mean square deviation have shown bigger values between the seepage model and the discharge hydrograph obtained with the nonlinear approach, as much as we varied the coefficient .
The distribution of flow between the matrix and fracture flow component of the system was estimated by applying the numerical model KarstMod. The results show that in the case of a linear model 19% of the total flow takes place in diffuse, and 81% in the fracture component of the flow. In the case of a nonlinear model, only 13% of seepage occurs in a diffuse phase, and 87% in the fracture phase. Variability in intensity of partial components of the flow in the linear variant is more pronounced than in the nonlinear as a direct consequence of a more intense variability of the input discharge hydrograph. Furthermore, the discharge hydrograph obtained by the KarstMod provides better results and more significant correlation with nonlinear model, compared to the linear one. This is probably due to greater hydraulic conductivity value of fracture () in the case of nonlinear approach which leads to more rapid response of aquifer to the rain. It should also be noted that the flow field, gained by the Modflow, shows velocities up to 0.26 m/s. Along the axis of the aquifer the average flow velocity in highly permeable fracture layer is 0.082 m/s, which is very close to the measured values of the apparent flow velocities in this area [16].
This paper shows that karst aquifer can be analyzed by using only limnigraphic and precipitation data and by taking appropriate assumptions in creating the conceptual model. Also, it is possible to identify some characteristic of the karst aquifer from the model simulation results, such as hydraulic conductivity coefficient. Furthermore, it is shown that this approach can be used in creating a proper transport model for the analyzed aquifer due to the relation between groundwater flow velocities and transport process, that is, advection component of transport, but greater investment in field research is certainly needed to define a model more similar to the conditions that exist in the natural environment. The application of a porous-equivalent media model Modflow to the Golubinka karst system indicated that for daily average hydrologic conditions the used method can meet calibration criteria for groundwater levels and average flows.
Conflicts of Interest
The authors declare that they have no conflicts of interest.