# Hydraulic-Hydrology Analysis of the Turbulent Seepage Flow within Karst Aquifer of the Golubinka Spring Catchment

Authors: Goran Lončar; Željko Šreng; Damir Bekić; Duška Kunštek
Journal: Geofluids
 (2018)

---

## 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.


---

## Body

## 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 (K), the hydraulic gradient (I), and the entire cross-sectional area of ​​flow (A), 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 P1; 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 (kt) 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 (Figure1). 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.Figure 1
Schematic hydrogeological map of the Golubinka spring catchment area.


## 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 datam 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 Q(t) for the Golubinka spring.Figure 2
Limnigraphic data on the water level difference between the “sea” and “freshwater” side of the partition wall built on the spring itself and the precipitation data from the omnigraphic station Zadar.To obtain the discharge hydrograph we introduced theKm coefficient of the Golubinka spring discharge and the water level difference m from the “sea” and “freshwater” side of the partition wall QGOLUBINKA(t)=f(m). The Km 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 QGOLUBINKA(t)=Km·m and the second assumes nonlinear correlation QGOLUBINKA(t)=Km·sqrt(m). The Km 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 (P1961–1990/P2012–2015=Q1961–1990/Q2012–2015). The calculated hydrographs for utilizing linear (for Km=0.4) and nonlinear (for Km=0.45) 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.Figure 3
Discharge hydrograph of the Golubinka spring in the analyzed period.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 P1 (Figure 1), located in the middle of the catchment, is the reference value for the entire catchment area.(1)Ht+Δt=Ht+-QGOLUBINKA-QPUMP+γ·P·RS·R·Δt,where Ht, Ht+Δt are groundwater level at the beginning and the end of the calculation [L], QGOLUBINKA is discharge at the Golubinka spring [L3T−1], QPUMP is pumping rate [L3T−1], γ is infiltration coefficient (adopted 0.6; [19]), R is catchment area (adopted 65 km2; [18]), and S 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 P1 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 P1 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 P1 (X = 15019′47,8′′E; Y = 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 P1. Figure 4 shows the calculated, (1), groundwater levels for the analyzed period 2012 to 2015 with adopted model parameters: R = 65 km2, S = 0.01, and γ = 0.3. The calculation of groundwater level obtained on the basis of discharge hydrograph with assumed linear correlation, QGOLUBINKA(t)=Km·m, yields the peak groundwater level of GWLmax(2012.–2015.)=61.29 m a.s.l. and the mean level of GWLmean(2012.–2015.)=48.11 m a.s.l. For the nonlinear correlation QGOLUBINKA(t)=Km·sqrt(m) the calculations yield higher peak level of GWLmax(2012–2015) = 62.35 m a.s.l. and higher mean levels of GWLmean(2012.–2015.)=49.04 m a.s.l. As previously stated, calculated values of GWLmax and GWLmean for both linear and nonlinear assumption are lower than the ones measured in 1967 (Table 1). The groundwater measurements on piezometer P1 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 R is divided by (i.e., for the number 2: γ=0.6; S=0.02; R=32.5 km2).Table 1
Difference between calculated (2012–2015) and measured (1967) values ofGWLmax and GWLmean.
ΔGWLmax=GWLmax(1967)-GWLmax(2012–2015) Δ G W L m e a n = G W L m e a n   ( 1967 ) - G W L m e a n   ( 2012 – 2015 ) QGOLUBINKA(t)=Km⋅m Q G O L U B I N K A ( t ) = K m · s q r t ( m ) Q G O L U B I N K A ( t ) = K m ⋅ m Q G O L U B I N K A ( t ) = K m · s q r t ( m ) 2.71 1.65 4.89 3.96Figure 4
Groundwater level for the analyzed period 2012–2015 and the position of the reference piezometerP1; calculations according to (1); red dot represents measured value (49.5 m a.s.l.).


## 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 formI∝vn. For the turbulent flow n=2 and the appropriate seepage law is v=ktsqrt(I), where kt 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: Q(y)=Ky·sqrt(I), where K(y) 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: K(L1,H1) and K(L2,H2). K was calculated using expression: K=kt·A, where A is cross section of the aquifer, A=H·B. By integrating the above-mentioned expression for turbulent flow in the fractures, and by separating the variables, (3) was obtained.(2)Kmean=KL1,H1+KL2,H22,(3)H1-H2+QQKmean2L=0,where Kmean is the arithmetic mean of the conductivity [L3T−1], K(L1,H1) and K(L2,H2) are hydraulic conductivities [L3T−1], L is distance between the edges of the section [L], H1 is groundwater level at piezometer P1 [L], H2 is groundwater level at Golubinka spring [L], and Q 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 kt (Figure 5), which represents the calibration parameter. The following values were used: ​​L = 10 km (distance P1 to Golubinka spring, Figure 1), BP1 = 2000 m (the width of the fracture system at the piezometer P1), BGolubinka = 150 m (width of the fracture system at the Golubinka spring), H1 is hP1(t) of Figure 4, and H2 is m(t) 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 kt-P1=1.0·10-6 m/s (at piezometer P1) and kt-GOLUBINKA=7.5·10-2 m/s (at Golubinka spring), and for the assumption of a nonlinear model the results were kt-P1=5.5·10-5 m/s and kt-GOLUBINKA=3.7·10-2 m/s.Figure 5
Time series of the amount of seepage through the Golubinka catchment (groundwater flow in the aquifer) and discharge at the Golubinka spring.
(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 (Figure6). 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 (QFS). The model has 3 water balance equations [23]:(4)dEdt=Peff-QEM-QEF,dMdt=QEM-QMS-QMF,dFdt=QEF-QFS,where E, M, and F are water levels in compartments E, M, and F, respectively, S is outlet of the system, Peff is the effective precipitation, assumed as 0.6·P [19], and QEM, QEF, QMS, QFS, and QMF are internal discharge rates per unit surface area.Figure 6
Structure of the model platform.Discharge at outletQS is defined as(5)QS=RQMS+QFS-QpumpS,where R is the recharge area of the catchment and QpumpS is the discharge rate abstracted at the outlet per unit surface area.Exchange functionQMF is defined as(6)QMF=kMF×sgn⁡M-F×M-FLrefαMF,where kMF is the specific discharge coefficient, αMF is a positive exponent, QMF is the algebraic flow from compartment M to compartment F, and Lref is a unit length.Input data of the model were precipitation, pumping rates from the Golubinka spring, and the discharge hydrograph from Figure3. 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 (BE).(7)NSE=1-∑US-Uobs2∑US-Uobs¯2,BE=1-∑QS-Qobs∑Qobs,where U stands for either the discharge or the square root of the discharge and Uobs 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 (0.5·NSE+0.5·BE), 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 criteriaWOBJ>(WOBJmin=0.6) are collected. Since more sets of parameters satisfied the criteria WOBJ>WOBJmin 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 WOBJ>WOBJmin 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 t, 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.Table 2
Ranges of parameter values carried out by the optimization procedure based on the input hydrograph obtained using (1) (linear and nonlinear) and the optimal values of the parameters used for calculating hydrographs with the KarstMod model.
Optimization parameter Range Optimal value Linear Nonlinear R (km2) 20–70 30.42 31.89 E 0 (mm) 0–100 1.76 3.66 M 0 (mm) 0–100 19.73 58.57 F 0 (mm) 0–100 25.98 18.92 k E F (mm/day) 1 E - 4–100 0.0551 1.37 α E F (-) 0.2–4 2.08 0.43 k E M (mm/day) 1 E - 4–100 0.0954 0.46 k M S (mm/day) 1 E - 4–100 0.0130 0.00035 k M F (mm/day) 1 E - 4–100 0.00026 0.00031 k F S (mm/day) 1 E - 4–100 0.0277 0.090 α F S (-) 0.2–4 1.00 0.68 WOBJ value (calibration stage) 0.61 0.69 WOBJ value (validation stage) 0.63 0.60Figure 7
Calculated discharges of the Golubinka spring obtained by the KarstMod model simulation, and the linear (Q=Km⋅m) and nonlinear (Q=Km⋅sqrt(m)) approach from the first step of the methodology.
(a) (b)Figure 8
Flow distribution between the diffuse (M – matrix) and the fracture (F) flow components through the analyzed karst aquifer.Figure 9
Q 0.05 S ( t ) and Q0.95S(t) discharge, 90% confidence interval for the simulated discharge (Q Golubinka–Karstmod); (a) linear, (b) nonlinear.
(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.(8)∂∂xkxx∂h∂x+∂∂ykyy∂h∂y+∂∂zkzz∂h∂z+W=SS∂h∂t,where kxx, kyy, and kzz are hydraulic conductivity coefficients in the x, y, and z coordinate axes, h is piezometric level, W is flow (volumetric) per unit of volume that represents the sinks and/or springs, the SS is coefficient of the specific storage in the porous media, and t 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 (Figure10) 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.Figure 10
The Modflow domain of the Golubinka basin and the location of Golubinka spring and piezometerP1.Parameterisation of the model in terms of determining appropriate values of hydraulic conductivity coefficientskx,y,z, the effective porosity neff, and total porosity coefficient ntot is adopted on the basis of the calibration procedures which rely on groundwater level at piezometers P1 (Figure 4). The adopted values are shown in Table 3 and Figure 11. The results for the groundwater level comparison for piezometer P1, 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 P1, during the simulation period 1.9.2012–05.26.2015.Table 3
Adopted hydraulic conductivity valueskx,y,z, the effective porosity coefficient neff, and the total porosity coefficient ntot, based on the calibration procedures relied on in the piezometer P1 groundwater level hydrograph.
Parameter Layer Surface Middle Bottom (variable) k x (m/s) 1.0 E - 4 1.0 E - 4 3.0 E - 1 / 1.0 e - 4 k y (m/s) 1.0 E - 4 1.0 E - 4 3.0 E - 1 / 1.0 e - 4 k z (m/s) 1.0 E - 4 1.0 E - 4 3.0 E - 1 / 1.0 e - 4 n e f f   ( 1 ) 0.02 0.02 0.01 n t o t   ( 1 ) 0.02 0.02 0.02Figure 11
Spatial distribution of hydraulic conductivity coefficientskx,y,z in the fracture layer with a thickness of 1 m (from −2 to −1 m a.s.l.). Blue area indicates the values kx,y,z=3.0E-1, and white kx,y,z=1.0E-4.Figure 12
Results of the piezometerP1 groundwater level hydrograph (applied coefficients kx,y,z, neff, ntot of Table 3).Figure 13
Hydraulic contour lines and the velocity vector in the bottom layer with a thickness of 1 m (from −2 to −1 m a.s.l.) on 15.11.2014, at the onset of the greatest groundwater level for piezometerP1 position 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(kt) 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 m 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 kt.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 (kFS) 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.



---

*Source: 6424702-2018-03-21.xml*