TIME SERIES ANALYSIS, MODELLING AND ASSESSMENT OF OPTIMAL exploitation OF THE NEMANJA KARST SPRINGS, SERBIA ANALIZA ČASOVNIH SERIJ, MODELIRANJE IN OCENA OPTIMALNEGA IZKORIŠČANJA KRAŠKIH IZVIROV NEMANJA V SRBIJI Igor JEMCOV^ & Metka PETRIČ2 Abstract UDC 556.34(497.11) Igor Jemcov & Metka Petrič: Time series analysis, modelling and assessment of optimal exploitation of the Nemanja karst springs, Serbia The time series analysis was applied in the case-study of a karst aquifer in Serbia in order to study its functioning, hydrody-namic behavior and hydraulic properties. Focusing on the definition of groundwater budget, due to very complex functioning of karst systems the correlation and spectral analyses were used to emphasize the importance of transforming the input data - precipitation to effective infiltration. The characterization of karst aquifer was further improved by separating the output component - discharge to base-flow and fast-flow components. Additionally, the importance of these transformations was proved in application of the regression model for the simulation of discharges based on the effective infiltration functions. A recharge-discharge model was applied in accordance with the active groundwater management, defining optimal "exploitable" regimes, which included the analyses of storage changes in karst water reservoirs under natural conditions, and calculation of the potential exploitation conditions. Keywords: karst hydrogeological system, time series, flow components, groundwater storage, exploitation capacity. Izvleček UDK 556.34(497.11) Igor Jemcov & Metka Petrič: Analiza časovnih serij, modeliranje in ocena optimalnega izkoriščanja kraških izvirov Nemanja v Srbiji Analizo časovnih serij sva uporabila pri proučevanju delovanja, hidrodinamičnih razmer in hidravličnih lastnosti kraškega vodonosnika v Srbiji. Pri oceni bilance podzemne vode sva zaradi posebnih značilnosti delovanja kraških sistemov z metodo korelacije in spektralne analize izpostavila pomen spremembe vhodnih podatkov iz padavin v efektivno infiltracijo. Karakterizacija kraškega vodonosnika je bila nadalje izboljšana z ločevanjem izhodne funkcije-pretoka v komponenti baznega in hitrega toka. Dodatno je pomen teh transformacij potrdila uporaba regresijskega modela za simulacijo pretokov na osnovi podatkov o efektivni infiltraciji. Model napajanje-praznjenje sva uporabila skladno z načeli aktivnega upravljanja s podzemnimi vodami in določila optimalne režime izkoriščanja. Pri tem sva analizirala spremembe uskladiščenja v kraških vodonosni-kih v naravnih pogojih in izračunala pogoje potencialnega izkoriščanja. Ključne besede: kraški hidrogeološki sistem, časovne serije, komponente toka, zaloge podzemne vode, kapaciteta izkoriščanja. INTRODUCTION Karst groundwater represents an important resource of water supply in many regions, Serbia being one of them. Water deficiency during the recession periods occurs as one of the main problems in water practice. Existing situation of water scarcity in Serbia may be significantly im- proved by means of adopting a proper procedure for the estimation of karst groundwater budget and storage. One of the frequently applied methods of rational management of the karst groundwater is to control its discharge regime by "borrowing" the water from the storage. The 1 University of Belgrade, Faculty of Mining and Geology, Department of Hydrogeology, Djusina 7, 11000 Belgrade, Serbia, e-mail: igor@jemcov.com 2 Karst Research Institute at ZRC SAZU, Titov trg 2, 6230 Postojna, Slovenia, e-mail: petric@zrc-sazu.si Received/Prejeto: 22.10.2009 main reason for opting for this method is that it enables a prompt compensation of water during the high water periods. Additionally, the "borrowing" is rational from the economic perspective. On the other hand, the major concern regarding this technique of tapping the karst groundwater rests with the possibility of overexploitation, which can further cause multiple problems (Pulido-Bosch 1999). Therefore, it is very important to establish the principles for karst aquifer characterization, and to properly assess the potential and available resources for groundwater tapping. The quantitative identification of karst hydrogeo-logical systems and the karst groundwater budget estimation are still underdeveloped fields of scientific study, mainly due to the complexity of the karst hydrogeologi-cal systems. Insufficient knowledge of the quantitative parameters is one of the principal problems when it comes to assessing the storage characteristics for the rational management of karst groundwater. Faced with the lack of data needed for the characterization of the water supply potential of karst aquifers, the analyses of spring hydrographs may provide valuable indirect information regarding the structure of karst hydrogeological systems (Mangin 1984). The assessment of the implementation potential of the measured regime control is an important issue and it serves as the basis for the preliminary studies. This is why a model which simulates the potential of the water exploitation has to be created and constantly improved in all phases of hydrogeological explorations. The paper intends to contribute to the study of functioning of karst aquifers, emphasizing the importance of transformation and separation of the recharge-discharge components in the groundwater budget equation. This was additionally confirmed in the application of the multiple regression model which simulates discharge rates based on the effective infiltration and established separated outflow components. Important objective of this paper is to develop the recharge-discharge model for simulating and assessing the optimal exploitation regimes of karst groundwater. GENERAL HYDROGEOLOGICAL CHARACTERISTICS OF THE STUDY AREA The Nemanja karst springs are located at Carpato-Balka-nides of Eastern Serbia (Fig. 1). Their catchment is mainly developed in Jurassic limestone and the karst aquifer shares the characteristics of a normal diffusive type of discharge. The main drainage points are the Nemanja and Klisura springs, all above the Mirosava riverbed level (Fig. 1). According Fig. 1: Simplified hydrogeological map of the Nemanja karst springs on a shaded relief model. legend: 1. Jurassic limestone-karst aquifer; 2. Tertiary sediments-intergranular porosity; 3. Permian sandstones-non-permeable rocks; 4. Alluvial depo-sits-intergranular porosity; 5. fault; 6. Captured karst spring; 7. Non-captured karst spring; 8. Ponor; 9. Main direction of karst groundwater flow; 10. verified connection between ponor and spring with dominant apparent flow velocity; 11. Surface flow; 12. Gauging station. to the system analysis approach, this is a typical binary karst system (Marsaud 1996) with a non-karstic part (Permian sandstones) providing allogenic recharge. The two karst springs, Nemanja 1 and Nemanja 2 (numbers marked within the symbols for captured and non-captured springs on Fig. 1), located at the contact with alluvial deposits and less permeable Tertiary sediments, represent low and high drainage points of a unique karst system. Their average annual discharge is 0.037 m3/s. During high waters discharge exceeds 0.22 m3/s and, at relatively fast pace, decreases to its minimal yield of 0.017 m3/s. The long-term mean annual precipitation is 600 mm. According to detailed hydrogeological exploration and hydrological budgeting, we assess the extent of the catchment to be somewhere below 6 km2. Its average altitude is less than 400 m asl. More than 30% of the karst area is covered by a thin layer of Tertiary sediments and vegetation. The vegetation index is 0.6. KARST GROUNDWATER BUDGET In order to rationally use and manage karst water resources, it is necessary to determine their overall water potential which is based on the existing groundwater budget. Due to a high complexity of the karst aquifer and data limitations, a system analysis approach was applied introducing a "karst hydrogeological system (KHS)". It represents a carbonate formation with developed karst porosity, which according to boundary conditions transfers and transforms the input components (fluids) to the output components. The input parameters are of a great importance, and they were determined by the means of correction of the obtained values of precipitation, including the interception (Petrič 2002). Measured precipitation data (P) were corrected (Pa) considering the errors due to aerodynamic effects and wetting loss (Sevruk 1982), and then reduced for the amount of precipitation intercepted by the vegetation cover (Int). To assess this amount the conceptual Rutter model was used (Rut- Fig. 2: Input parameters (measured and transformed) in KHS, and hydrographs of discharge and separated components of the Nemanja karst springs. Legend: P. Measured precipitation (a); Pa-int. Corrected precipitation, reduced for the amount of precipitation intercepted by the vegetation cover (b); SWE. Snow-water equivalent (c); Pa-int-SWE. Actual precipitation that reaches the ground (d); Ref. Effective infiltration (e); Qsum. Summary hydrograph of the Nemanja 1 and Nemanja 2 karst springs; Q N1. hydrograph of the Nemanja 1 karst spring; base-flow. hydrograph of the base-flow component of the summary outflow of Nemanja 1 and Nemanja 2 springs (f). ter et al. 1971). The relevant analysis of the snow melting processes was especially taken into account (US Army Corps of Engineers 1998) and the snow-water equivalent (SWE) was assessed. The effective infiltration (Ref) of water into the karst hydrogeological system, as well as the evapotranspiration (ETR), were ascertained through the analysis of the soil budgeting (Reed 2003). An important characteristic of karst aquifers is duality of groundwater flow (White 1969; Atkinson 1977) and based on this principle most conceptual models consider two types of flow, e.g., base-flow and fast-flow. Study of the relationship between these two types of in- terconnected flows provides valuable information about the hydrodynamic behaviour of karst aquifers (Jem-cov 2007). For this purpose, a computerized base-flow method of Institute of Hydrology (1980a, b) was used for the karst hydrograph separation (Fig. 2). For the turning point the test factor f = 0.9 was used, while the parameter N (the number of days over which a minimum flow is determined) was assessed according to the cross-correlation function of fast-flow with a relatively fast drop below the level of significance (Jemcov & Petric 2009). For the Nemanja springs, the share of the base-flow in the total flow is 87%. TIME SERIES ANALYSIS Karst aquifers are characterized by high heterogeneity and spatial variability of hydrogeological parameters. The hydrograph analysis is often applied with a view to study the functioning and hydrodynamic behavior of the karst aquifer (Bonacci 1993; Kresic 1997). Valuable indirect information regarding karst systems can be obtained as a result of the time series analysis (Box & Jenkins 1970; Yevjevich 1972). Mangin (1984) developed a special methodology of study of the input-output relations in the karst aquifers. This methodology was based on the systemic approach, and it was applied and further developed by a number of authors (Padilla & Pulido-Bosch 1995; Larocque et al. 1998; Panagopoulos & Lambrakis 2006; Jemcov & Petric 2009; etc.). The autocorrelation of the flow rates of Nemanja springs (Fig. 3a) exceeds the confidence limits for approximately 122 days, which implies that storage is significant and water is released from the aquifer gradually. Additionally, the slope of the autocorrelogram initially drops quickly (for less than 10 days), and afterwards more slowly. This bimodal behavior indicates the duality of the karst aquifer (Panagopoulos & Lambrakis 2006). The autocorrelogram of the fast-flow component of the Nemanja springs drops below the level of significance in 18 days, which confirms previous statement and indicates a noticeable influence of the fast-flow component on the total outflow. The autocorrelogram of the base-flow component of the Nemanja springs diminishes very gradually and has a significant influence on the total outflow. The spectral density function of daily discharges of the Nemanja karst springs (Fig. 3b) shows high peaks at the low frequency of 0.0027, which confirms the presence of important annual cycle. Additionally, there are several peaks of low densities up to the frequency of 0.05, and at the frequencies higher than 0.15 the function inclines to zero. The regulation time Treg is 73 days, which indicates a very long impulse response, which implies a significant storage and well structured KHS (Larocque et al. 1998). Considering the spectral density function of the base-flow, this component has a significant influence on the summary outflow at low frequencies, which corresponds to the peaks of 41 days. Afterwards, at the frequencies which correspond to the period of 38 days, the fast-flow component takes a complete control of the spring discharge. The cross-correlation function (CCF) of the Nemanja karst springs (Fig. 4a) shows non-symmetrical behavior and low level of influence of precipitation on discharge rate. This function becomes insignificant after 12 days, and afterwards (64 days) it exceeds the level of significance as the consequence of a random behavior of the input component (precipitation). A considerable increase of the correlation coefficient was achieved by the change of the input variable from precipitation to the effective infiltration. The increase in the values of r(k) is gradual - it is the highest for the effective infiltration, which becomes insignificant after 137 days. Obtained results confirm the significance of the transformation of the input components. Similar improvements were achieved when it comes to the base-flow component. When we compare the effective infiltration and the base-flow (Fig. 4b), the CCF slowly decreases with the peaks at 7, 50, 80 and 120 days, which implies to non-homogeneous KHS, e.g., to different responses in various parts of the system or the existence of annexes to drainage systems - ADS (Mangin 1975). A significant attenuation of the impulse response of the base-flow CCF is highly controlled by the characteristics of the karst hydrogeological system, fig. 3: Autocorrelograms (a) and spectral density functions (b) of the Nemanja KHS. legend: Qn. Summary discharge of the springs; Qb. Base-flow component; Qf. fast-flow component; Ref. Effective infiltration; ls%. level of significance. Note: for (Qn) the spectral density function exceeds 60,000, and for (Qb) it exceeds 40,000; these values are not shown on the graph. which is the consequence of a highly structured part of KHS and indicates slow travel of the pressure pulse through the top soil and unsaturated zone, with significant lateral component in the epikarst zone. Additionally, influence of the snow-melting on the base-flow component is noticeable. The CCF for the fast-flow component shows synchronized behavior with the CCF of the summary outflow, with a relatively high level of dependence of the input-output components (Fig. 4c). Relatively significant differences in the CCF of input parameters ((P-Int-SWE) and Ref) and base-flow com- ponent are the consequence of an influence of the soil moisture, which serves as an important filter and transforms the input parameters. The same influence is reasonably lower for the fast-flow component, since this part of the infiltrated water was not strictly connected to the condition of soil moisture capacity. Namely, some previous detailed researches of the infiltration process have demonstrated that the so-called rapid infiltration may take place entirely independent of the soil moisture condition when the precipitation water infiltrates through the cracks in the soil (Rushton & Ward 1979; Fig. 4: Cross-correlograms of the Nemanja springs. legend: Qn. Summary discharge of the springs; Qb. Base-flow component; Qf. Fast-flow component; P. Measured precipitation; P-Int. Corrected precipitation, reduced by the amount of precipitation intercepted by the vegetation cover; P-Int-SWE. Actual precipitation that reaches the ground; Ref. Effective infiltration; Ls%. level of significance. fig. 5: Plots of functions of cross-amplitude (a); coherency (b); obtained for the studied KhS. Petrič 2002). Obtained results confirm the previous conclusions about a well-structured karst hydrogeological system (KHS) with prevailing shaft flow in the epikarst zone (Klimchouk 2004), and transmission of a pressure pulse trough the saturated zone functioning as a water hammer (Larocque et al. 1998). ^e cross-amplitude function (CAF) shows (Fig. 5a) that KHS notably filters and attenuates the input signal at high frequencies and increases it at low frequencies. The observed peak at low frequencies clearly indicates periodicity, e.g., annual and seasonal cycles. tte coherency function (COF) shows (Fig. 5b) a non-linear relation between the input and output variables with an average value of approximately 0.6. tte COF for Ref reveals a linear character at low frequencies (0-0.003), with oscillatory values in the range of 0.43-0.99 (mean 0.82); while the average value for the whole range is 0.69. Fig. 5: Plots of functions of gain (c); and phase (d); obtained for the studied KHS. The gain function (GAF) expresses (Fig. 5c) the relation between the base-flow and the fast-flow (Padilla & Pulido-Bosch 1995). Dependence between (Pa-Int-SWE) and Qn is characterized by an intensive alteration of the base-flow and fast-flow, forming an intermediary flow which exists in a wide frequency range (0.05-0.5). Results of the GAF for Ref-Qn relationship are significantly different and confirm previously explained influence of the base flow in KHS of the Nemanja springs. The obtained values of the phase function (PhF) were mainly non-coherent and unsorted in a wide range of frequency (Fig. 5d), except when it comes to Ref as the input data, where the delay of approximately 45 days was registered. Observed linear character between the input (Ref) and output (Qn) components at low frequencies (< 0.003) in the COF (Fig. 5b), GAF (Fig. 5c) and PhF (Fig. 5d) presents a rapid response of KHS (i.e., heavy storms, with prevailing shaft flow). Out of this frequency range (i.e., regular rainfalls), the KHS responses to the inputs in a non-linear character. SIMULATION MODEL Obtained characterizations of the KHS confirm the importance of transformation of the input components (recharge) for the two types of interconnected flows (base and fast). According to this relation, the recharge-discharge model is formed. ttis model is not strictly stochastic, since it takes into account the physical processes and the nature of the system. ttis is why it should be understood as a stochastic-conceptual model (Kresic 1997). Using this model, the input component (Ref) is transformed into the components which simulate the functioning of KHS; introducing two functions - Kfun-base and Kfun-fast - which are related to the base-flow and fast-flow components. Schematically observed, the outflow is not depending only on the infiltration values - it also depends on previously accumulated water (storage) in the KHS. For that reason, a fictive state of storage - vs0_i is simulated, as a cumulated value of the previous state of storage. tte efi'ective infiltration - R0_i is reduced for the value of fictive outflow of the previous day - Kfun-base/fast0_i. tte value of fictive outflow is a result of the values of fictive state of storage, multiplied by k coefficient, which simulates a fictive depletion correlated with the base-flow (kb) and the fast-flow (kf). tte fictive depletion coefficients are estimated by the calibration process with a view to reach a maximal correlation coefficient with the base-flow and fast-flow. tte initial value of the state of storage is also obtained during the process of calibration (Tab. 1). To simplify the parameter estimation of the state of storage, the periods of hydrological years (i.e., end of recession period) are analyzed and in this way the estimation of the initial values of storage for the fast-flow addition, during the process of calibration of both parameters (storage and fictive depletion), the estimated values of groundwater budget were also considered in the model. The established relation of base-flow components and Kfun-base parameters shows relatively high dependence, while dependence of fast-flow components and Kfun-fast parameters imply a wide dispersion of data (Fig. 6). Considering a non-linear character of the input-output relationship of KHS, multiple regression was applied for the purposes of the discharge simulation Q=f(Kfun-base, Kfun-fast) (Fig. 7). Established function of regression is polynomial and it is shown below: y = a + b • + c • x2 + d • x3 + e • x4 + f • x2 + ... ... g • xi + h • x3 + i • x4 + j • x5 Where x1=Kfun-base-; x2=Kfun-fast, y=Q. With function parameters a b c d e 1.358 -7.153 1.886 -5.822 1.036 f g h 0.705 -0.417 0.218 i J -3.4x10-2 1.6x10-3 Considering the fact that what was applied is a simple model, one can state that a good correlation coefficient (r-0.84) between measured and simulated val- Tab. 1: Scheme of calculation of the fictive outflow Kfun-base/fast. R1 R2 R3 Vs0 = ? Vsi = Vs0 + Refi - Kfun0 Vs2 = Vsi + Ref2 - Kfuni Vs3 = Vs2 + Ref3 - Kfun2 Kfun - base / fast0 = Vs0 x Xb/f Kfun - base / fasti = Vsi x Xb/f Kfun - base / fast2 = Vs2 x Xb/f Kfun - base / fast3 = Vs3 x Xb/f Qbase/fast0 Qbase/fast1 Qbase/fast2 Qbase/fast3 Ri —► Vsi = Vsi-1 + Refi - Kfuni-1 Kfun - base / fasti = Vsi x Xb/f R - correlation coef. Kfuni - Qbase/fasti = max Qbase/fasti component is avoided. Furthermore, the calibration of this parameter was achieved considering similar values at the start and at the end of the observed period. In ues was obtained. We can observe a high level of correspondence at the recession part of the hydrograph, and a low level of correspondence at the peaks of the hydro- 1J5 14 a 1.05 0.35- OZ 0.4 0.6 08 Qbase i ■y ■•y . ■■ ■ •■■■■■ i " ' y \ /i Pi . 1/ - /A j''*- -i.. C 12 fig. 6: Dependence of the flow components (base and fast) and the estimated Kfun-base and Kfun-fast parameters, obtained as a result of the effective infiltration of Ref multiplied by the fictive coefficient of depletion Xb and Xf (dotted line-confidence limit of 95%). Fig. 7: 3D graphical representation of applied discharge function of the Nemanja springs. (x1=Kfun-base, x.,=Kfun-fast, y=Q, cross marks-input data, surface colour area-established model). graph (Fig. 8). The latter is caused by linear and instant responses of the system to heavy storms (Figs. 4a, 4c, 5b, 5c and 5d). Furthermore, it could be a consequence of a non-established relation for the epikarstic zone (i.e., simulating piston effect), and many uncertainties related to the turbulent flow. fig. 8: hydrograph of simulated (Qsim) and measured discharges (Qn) of the Nemanja karst springs. OPTIMAL EXPLOITATION CAPACITY Quantitative analysis of a karst hydrogeological system facilitates a detailed determination of the ways in which the karst systems function, which makes it easier to choose the right technical solution of the tapping structure. Characterization of the Nemanja karst springs shows a good potentiality of groundwater management, based on the principle of "water-borrowing" from the storage in the times of recession. The main reason for selecting this concept is relatively fast replenishment of water during the high-water periods, which is normally accomplished in one hydrological cycle. If we want to apply this concept, it is of the utmost importance to determine the optimal exploitation capacity in the first place. The optimal capacity defines a sustainable water exploitation amount under different conditions, respecting water demands, ecological criteria etc., and avoiding overexploitation. Prior to defining the "exploitation" regime, a change of storage in KHS must be completed under natural conditions. Based on the inflow and outflow relations, the values of changes in storage (AV,) in the studied KHS were collected by adding and varying the continuous cumulative values to the initial storage (Fig. 9). ttus, information on the state of water storage (karst accumulation) in the natural conditions is gathered via this relation: Vi=AVi+Vi_1. tte most sensitive part of this analysis is the estimation of initial storage. Knowledge of this particular parameter is strictly connected to the level of hydrogeological exploration, which should not be neglected, and could lead to underestimation or overestimation of possible exploitation capacity (Jemcov 2007). Two different scenarios of storage changes can be expressed through the analysis of potential exploitation: Qpi > Qei Vpi = Vei AQe'i = 0 AVe'i = 0 Qp Vpi (1) i = 1, 2.....n - days Qp - discharge of karst spring under natural conditions; Vp; - state of storage in natural (unchanged) conditions of discharge regime; Qe; - simulated discharge regime under condition of exploitation; Qpi - fictive (changed) discharge which will appear at karst spring as the consequence of changes in the state of storage in case the exploitation is instantly canceled; Qei - fictive discharge of karst spring; Vei - state of storage in KHS under the exploitation regime, which corresponds to fictive discharge - Qp,' ; I - perod of "water-borrowing" from the storage; t - period of cumulated water in storage. fig. 9: Estimated state of storage (a) and simulated exploitation conditions (b) in KHS of the Nemanja springs. Fictive discharge of karst spring corresponds to natural (unchanged) spring discharge (Q ), only if AQei=0 and AVe;=0, ^ and otherwise in case: Qpi < Qei Vpi < Vei Qe'i Ve'i i (2) When the process of "water-borrowing" from the storage in KHS is started, it provokes cumulative lowering of the storage state. On the contrary, the newly-infil- trated water increases the state of storage until it reaches the condition referred to in the Equation (1). The state of storage in KHS in the condition of exploitation can be portrayed by the following relation: Vei = Vpi - (Qei - Qpi) + iQPi-1 - Ve'-i) curent_wat^■^_defic;t ^'rev;ously_ fo'rmed_water_defic;t (3) According to the previous relation (Eq. 3), the state of storage in KHS in condition of exploitation is actually the difference between natural conditions on the one hand, and cumulated "borrowed" water on the other hand. Here, the main obstacle may be the fact that it is difficult to determine the fictive discharge of karst spring - Qp'. tte changed (exploitation) condition of the state of storage will lead to the change - or, more precisely, to the lowering - of the discharge. This can be resolved by introducing the (fictive) coefficient of recession (a'), which is in inverse proportion to state of storage (Milanovic 1981): qp; Vei For the purposes of defining the fictive discharge, we assumed that KHS will keep the same recession conditions as in its natural (unchanged) state. Under this presumption, the values of fictive recession coefficient can be obtained as mean value of every recession episode. Using the estimated water resources stored in the karst aquifer under natural conditions as a baseline point, further analyses were conducted to assess variations in storage under artificial conditions (exploitation potential, limits and optimal values). According to adopted coefficient of seasonal variation of exploitation (1.4), the optimal capacity was reached at these values: 17.6-30.8 l/s (Fig. 9). By means of the aforesaid model, we can define the long-term optimal variant of exploitation. This model does not examine the type of artificial karst aquifer regulation, but simply defines quantitative conditions of possible rational exploitation, with the aim of satisfying the needs of the users in a multi-year period. (4) CONCLUSIONS Quantitative analysis and budget analysis of the elements of karst hydrogeological system (KHS) are both complementary to the traditional techniques of the hydro-geological exploration. They represent a valuable tool of integrated management of the karst aquifer. Transformation of precipitation data into effective infiltration emphasizes the importance of the structure of karst system, and prevents the influence of other processes (such as meteorological parameters, vegetation or soil). Our ability to determine the relation between the two types of interconnected flows (base-flow and fast-flow), provides us with the valuable information on the storage capacities of the aquifer, which is the essential condition for effective groundwater management. tte importance of duality of groundwater flow was confirmed by the use of the simulated model. The applied concepts that are based on karst groundwater budget provide valuable information on storage changes in KHS; they enable future predictions regarding the optimal exploitation rates, and facilitate karst water management. The results that are obtained by this analysis should contribute to the feasibility studies, and help us avoid the problem of overexploitation. ACKNOWLEDGEMENTS The authors are grateful for the valuable review comments and suggestions from the editor Nico Goldschei-der and three anonymous reviewers. REFERENCES Atkinson, T.C., 1977: Diffuse flow and conduit flow in limestone terrain in the Mendip Hills, Somerset (Great Britain).- Journal of Hydrology, 35, 93-103. Bonacci, O., 1993: Karst springs hydrographs as indicators of karst aquifers.- J. Hydrol. Sciences, 38, 5162. Box, G.E.P. & G.M. Jenkins, 1970: Time series analysis: Forecasting and control.- Holden-Day, pp. 575, San Francisco. Institute of Hydrology, 1980a: low Flow Studies.- Wallingford, Oxon, United Kingdom, Report No. 1, 1-41. Institute of Hydrology, 1980b: low Flow Studies.- Wallingford, Oxon, United Kingdom, Report No. 3, 12-19. Jemcov, I., 2007: Water supply potential and optimal exploitation capacity of karst aquifer systems.- Environmental Geology, 51, 5, 767-773. Jemcov, I. & M. Petric, 2009: Measured precipitation vs. effective infiltration and their influence on the assessment of karst systems based on results of the time series analysis.- Journal of Hydrology, 379, 304-314. Klimchouk, A., 2004: Towards defining, delimiting and classifying epikarst: Its origin, processes and variants of geomorphic evolution.- Re-published (modified) from: Jones, W.K. et al. (eds.). Epikarst, Proceedings of the symposium, 1s'-4'h October, 2003, Shepherd-stown, West Virginia, USA. Karst Water Institute Special publication, 9, 23-35, Charles Town. Also available at: http://www.speleogenesis.net. Kresic, N., 1997: Quantitative solutions in hydrogeology and groundwater modeling.- CRC Press LLC, pp. 480, Boca Raton. Larocque, M., Mangin, A., Razack, M. & O. Banton, 1998: Contribution of correlation and spectral analyses to the regional study of a large karst aquifer (Charente, France).- Journal of Hydrology, 205, 217-231. Mangin, A., 1975: Contribution ä l'etude hydrodynamique des aquiferes karstiques.- ttese Univ. Dijon. Annales de speleologie, 29/3: 283-332, 29/4: 495-601, 30/1: 21-124. Mangin, A., 1984: Pour une meilleure connaissance des systemes hydrologiques a partir des analyses cor-relatoire et spectrales.- Journal of Hydrology, 67, 25-43. Marsaud, B., 1996: Structure et fonctionnment de la zone noyee des karst ä partir des resultats experimentaux.-PhD ttesis. Paris XI, Orsay, pp. 305. Milanovic, P., 1981: Karst Hydrogeology.- Water Resources Publications, pp. 434, Littleton. Padilla, A. & A. Pulido-Bosch, 1995: Study hydrographs of karstic aquifers by means of correlation and cross-spectral analysis.- Journal of Hydrology, 168, 73-89. Panagopoulos, G. & N. Lambrakis, 2006: The contribution of time series analysis to the study of the hy-drodynamic characteristics of the karst systems: Application on two typical karst aquifers of Greece (Trifila, Almyros Crete).- Journal of Hydrology, 329, 698-376. Petrič, M., 2002: Characteristics of recharge-discharge relations in karst aquifer.- ZRC Publishing, pp. 154, Postojna-Ljubljana. Pulido-Bosch, A., 1999: Karst water exploitation.- In: Drew, D. & H. Hötzl (eds.) Karst hydrogeology and human activities. Impacts, consequences and implications. International Contributions to Hydrogeol-ogy, Vol. 20. Balkema, pp. 225-234, Rotterdam. Reed, S.M., 2003: Soil-Water Budget. Part I. Methodology and results.- FAO/UNESCO Water Balance of Africa Project. [Online] Available from: http:// www.ce.utexas.edu [Accessed 18*^ May 2009]. Rushton, K.R. & C. Ward, 1979: tte estimation of groundwater recharge.- Journal of Hydrology, 41, 345-361. Rutter, A.J., Kershaw, K.A., Robins, P.C. & A.J. Morton, 1971: A predictive model of rainfall interception on forests. I. Derivation of the model from observations in a plantation of Corsican pine.- Agricul. Meteorol. P, 367-380. Sevruk, B., 1982: Methods of correction for systematic error in point precipitation measurement for operational use.- WMO, No 589, Operational Hydrology, Report No. 21, Geneva, 1-91. U.S. Army Corps of Engineers, 1998: Engineering and Design - Runoff from Snowmelt.- Department of the Army, Washington, CECW-EH Engineer Manual 1110-2-1406. White, W.B., 1969: Conceptual models for carbonate aquifers.- Ground Water, 7, 3, 15-21. Yevjevich, V., 1972: Stochastic Processes in Hydrology.-Water Resources Publications, pp. 302, Fort Collins.