Definition of the river Gacka springs subcatchment areas on the basis of hydrogeological parameters

The river Gacka springs catchment area is located in the Dinaric karst, which is globally known as the locus typicus, or classical Karst. It is composed of four major and several minor karst springs of different discharge rates. The river Gacka springs are characterised by great discharge and exceptional quality, so the catchment area of the river is indicated in the Water Management Strategy (OFFICIAL GAZETTE NO. 91/08) as an area with strategically important reserves of drinking water for the Republic of Croatia. To determine the hydrogeological characteristics of the subcatchments of this large and complex aquifer system, hydrological and hydrochemical parameters were measured on the main springs. Data collected on the springs were analysed using the recession analysis by the „matching strip“ method, the statistical analysis of a time series of measured data both by autocorrelation (analysis of individual series) and by cross-correlation methods (analysis of interrelationships between time series), multivariate statistical analysis (Factor Analysis) of hydrochemical parameters using the software package STATISTICA 6.0 (1998), and geochemical modelling of hydrochemical parameters using the NETPATH computer program. Interpretation of lithological, structural and tectonic characteristics of the rocks, together with tracing data and the applied analytical methods, allowed the springs catchment of the river Gacka to be divided into three subcatchments. The results of this study imply the necessity of a multidisciplinary approach to research.

The catchment of the river Gacka springs is located in the Croatian part of the Dinaric karst, which represents the karst locus typicus.The catchment consists of four major (Tonković, Majerovo Vrilo, Klanac and Pećina) and several minor springs (Jaz, Marusino Vrilo and Graba) of different discharge rates.The Gacka springs are characterised by a high discharge and exceptional water quality, which was why this catchment area was proclaimed by the Water Management Strategy (OFFICIAL GAZETTE NO. 91/08) as one of the strategically important area of drinking water reserves of the Republic of Croatia.The research area is characterised by features of typical karst geomorphology.Research and data collection is difficult due to the heterogeneity of the karst aquifer, as well as the unknown geometry of the voids through which the water flows underground to the springs.
In addition, with a size of approximately 500 km 2 the river Gacka springs catchment area is categorised as a "big catchment" (KENDALL & MCDONNELL, 1998), which further complicates the study.Results of previous hydrological and hydrogeological studies (PAVIČIĆ et al., 2003;LUKAČ RE-BERSKI, 2008) have shown that the river Gacka springs catchment area can be divided into subcatchments of the main springs.Former detailed hydrological analyses were performed on the basis of hydrological data measured only on the profiles of the river Gacka (Bonacci & Andrić, 2008;Lukač Reberski, 2008).Hydrograph analyses performed based on these data are related to the reaction of the entire aquifer system to precipitation, while the characteristics of particular sub-parts of the catchment remain unknown.Therefore, detailed hydrological analyses were performed on the most important springs (Tonković, Majerovo Vrilo, Klanac and Pećina), providing new insight on this complex and strategically important karstic hydrogeological system.The investigations and the results presented here are part of an extensive and systematic hydrogeological study described within the authors' PhD thesis (LUKAČ REBERSKI, 2011).

GEOGrAPHIC, CLIMATOLOGIC AND GEOMOrPHOLOGIC CHArACTErIsTICs
The investigated area can be described as a hilly karst area with a large karst polje (Gacka polje) surrounded by the Ve lebit, Kapela, and Plješivica mountains.The highest altitudes (<1200 m.a.s.l) belong to the central part of the karst hinterland, whereas the river Gacka springs formed at the lowest altitude, (~450 m.a.s.l).Locally, the relief plays a key role in the definition of climate conditions.High mean values of annual air temperature are the highest in Croatia.There is also a significant difference in the spatial distribution of precipitation.These, are the direct result of the orographic characteristics of the study area.The mean annual air temperature ranges between 4 and 9°C, (ZANINOVIĆ et al., 2004) while the annual mean precipitation values range between 1000 and 2500 mm (GAJIĆČAPKA et al., 2003).The research area is characterised by typical karst geomorphology i.e. karst features are distinguishable in the relief (karrens, sinkholes, caves, karst poljes, karst springs, estavelles, sinks).

GEOLOGICAL, HYDrOGEOLOGICAL AND sTrUCTUrAL-TECTONIC CHArACTErIsTICs
The study area is composed predominantly of carbonate rocks with fractured and fracture-cavernous porosity.Over the wider area of the river Gacka springs catchment, various Jurassic and Cretaceous limestones and dolomites predominate.The great permeability of the limestone is the direct result of its fragmentation caused by intense tectonics, as well as by the lithological constitution which enables dissolution processes.Besides the carbonates, Tertiary clastites, i.e.Jelar deposits (BAHUN, 1974) are also abundant in this area, and can be observed overlying the carbonates, sometimes as only a few metre thick erosional remains, or up to several hundred metre thick deposits.These deposits are represented by limestone breccias made of unsorted Jurassic, Cretaceous and Palaeogene sediment fragments (VLAHO VIĆ et al., 2009).The Jelar formation has specific hydrogeological characteristics and, due to its lithology and spatial distribution, a specific hydrogeological function.The presence of marly-matrix breccias intercalated with marl lenses led to two opposing hydrogeological effects i.e. firstly the unusually well-developed surface and underground karst phenomena, and secondly the formation of relatively less permeable environments (BAHUN & FRITZ, 1975).These deposits can be observed at the bottom of the polje near the river Gacka spring, as well as in the western part of the catch ment.Fine-grained sediments with intergranular porosity occur in the form of the Quaternary cover in poljes and depressions and they do not have an important hydrogeological influence.
Lithology had an essential role in the formation of the hydrogeological properties, especially permeability (LUKAČ RE BERSKI et al., 2009).Hence, based on their permeability the deposits are divided into: -Very permeable carbonate rocks, which are the most common types, and where limestones of Jurassic and Cretaceous ages prevail; -Predominantly permeable rocks are represented by other Jurassic and Cretaceous carbonates and Palaeogene clastites, i.e. the Jelar deposits.The Jelar deposits exhibit about 50% less permeability than the Jurassic and Cretaceous carbonates, which surround them (BAHUN & FRITZ, 1975); -Predominantly impermeable rocks -comprising of dolomites of Jurassic and Cretaceous periods; -Impermeable deposits include the Eocene marls occurring in a very small region in the eastern part of the investigated area; -Variably permeable deposits are Quaternary deposits with heterogeneous properties, as their permeability varies depending on the thickness and composition.They are represented mostly in the depressions of karst polje.
One of the most important, and currently unsolved, karst problems is determining the boundaries and the surface of the catchment, for which classical geological and hydrogeological approaches are necessary but not always sufficient.Hydrological approaches may provide answers to some questions regarding the size of the catchment surface but cannot imply the position of the boundaries (BONACCI, 1995).Unlike a topographic water divide, determination of the groundwater divide in karst terrain depends on many factors.Furthermore, the groundwater divide is not constant as its po sition changes depending on the groundwater level (ŽUGAJ, 2000).Therefore a multidisciplinary research approach and data collection is required in order to precisely determine the location of the divide and the size of the catchment.Several tracings tests were carried out using the Nafluorescein artificial tracer to determine the apparent velocities and groundwater flow directions as well as the size of the catchment area.Tracing test results, performed between 1957 -2010., were presented in unpublished technical reports.Recent investigations suggest that the surface of the Gacka catchment up to the Čovići profile extends to 516 km 2 (BIONDIĆ et al.

2010).
Groundwaterflow directions depend more on the structural relationships than on the lithological characteristics because of the relatively monotonous lithology.The study area is characterized by the NW-SE Dinaric strike of the structures.Tracings showed that the groundwater flow direction in the greater part of the catchment area is parallel to these structures.This is due to the position of the main boundary faults, favourable to the stress direction (20°-45°), and along which the right transcurrent shifts appear, apart from the opening up of spaces and appearance of structures of a pullapart type (PRELOGOVIĆ, 1989).A fault direction unfavourable to the stress direction leads to local compression by closing of the space and possible prevention of water flow (Fig. 1).

HYDrOLOGICAL CHArACTErIsTICs
The river Gacka is characterized by favourable hydrological conditions.The flow regime is standardized compared to other karst rivers, and considerable variations of the flow quantities through its river bed do not occur.The ratio of the lowest, medium and highest discharge is 1:4:20, whereas the ratio of the neighbouring Lika river is 1:100:800 (BONACCI, 1987), which shows its vividly torrential character.
During the monitoring period (2008 -2010) the amount of precipitation varied, with 2008/09 being considered as an average year when compared to the long term average, where as 2009/10 was wetter.The average annual discharge during the monitoring period, (determined from data collected on the Gacka river gauge at Čovići Podgora, Fig. 1), was 13.44 m 3 s -1 of which 29% is related to the Tonković spring, 26% to the Klanac spring, 24% to the Majerovo Vrilo spring, 15% to the Pećina spring and 6% to the remaining springs (Fig. 2).
Recent investigations of the hydrological balance of the catchment indicated the significant retention properties of the karst underground system in the river Gacka catchment .1-very permeable rocks, 2-predominantly permeable rocks, 3-predominantly impermeable rocks, 4-most important faults, 5-extension zone, 6-swallow hole (ponor) and vertical cave (pit), 7-spring: permanent and intermittent, 8-groundwater divide, 9borehole, 10-groundwater connection (tracer test), 11-river, 12-river and rain gauge stations, 13-groundwater velocity determined by tracing tests (cm/s).(LUKAČ REBERSKI, 2008).Therefore the water balance of the catchment should be prepared over longer time periods because it provides more reliable discharge coefficient information.Due to the inadequate distribution of rain gauge stations (Fig. 1) in the river Gacka springs catchment, mean an nual precipitation was calculated using the regression equation defined for the research area (GAJIĆČAPKA et al., 2003).This equation also takes into account the precipitation change with altitude as all the rain gauge stations in Lika are below 760 m.a.s.l. and most of the catchment is located at higher altitudes.Mean annual precipitation is 1383 mm while the mean catchment discharge is determined from 30 years of daily discharge data (1972 -2002), so the mean annual discharge in the profile of the hydrologic gauge station ČovićiPodgora is 14.2 m 3 /s.
The discharge coefficient c can be defined as follows: (1 where Q is the mean discharge (m 3 s -1 ), T is duration of the mean discharge (s), P is the amount of precipitation (m), and A is the surface of the catchment (m 2 ).
The discharge coefficient of the catchments calculated using the above formula is 0.63, which means that 63% of total precipitation in the catchment infiltrates into the ground and flows towards the river Gacka springs.

METHODs AND rEsEArCH TECHNIQUEs TO DETErMINE THE DYNAMIC CHArACTErIsTICs OF THE sPrINGs
To determine the dynamic characteristics of individual springs and their catchments, daily water discharge rates were collected between 2005 -2010 from the Meteorological and hydrological service of Croatia (DHMZ).Due to the extreme ly heterogeneous characteristics of the karst aquifer, spring runoff hydrographs were used because they refer to the total response of the aquifer to precipitation (PADILLA et al., 1994).Analysis of spring hydrographs is both a cheap and very useful method in hydrogeological research, especially for karst aquifers.
Both the recession parts as well as the complete hydrographs were analysed.Recession analysis refers to the analysis of the falling limbs of the runoff hydrographs.The advantage of the recession method is that it does not need prior knowledge of the distribution of the potential and the individual parameters of the aquifer, although the recession coefficient directly depends on them.Recession curves mostly represent only part of the dynamic reserves of the aquifer, depending on the initial water level.The total water volume stored inside the aquifer can be much higher.Numerous authors have publicised the usefulness of this method (TAL-LAKSEN, 1995;DEWANDEL et al., 2003;BAKALOW-ICZ, 2005).In the case of long recession periods, which a lack of heavy and long precipitation periods, the recession curve could provide good insight into the structural characteristics of the aquifer.Both the shape and the characteristics of the recession curve, i.e. the calculated recession coefficient (α) values, imply the hydrogeological characteristics of the aquifer and are dependent on various factors, among which the most important are: porosity type, hydrological conditions, i.e. the water level as well as the rate of water supply from other aquifers.If the α values are larger than 10 -2 a rapid drainage system is indicated, with large fractures and channels; when the α value is smaller than 10 -2 , slow drainage from smaller pores and fractures may occur, i.e. it indicates that the rock matrix has the dominant role in the water flow through the karst underground (KREŠIĆ, 1997).
For recession analysis, it is of major importance that the recession period lasts for a long time period, at least for several months.Such a case has never been registered for the region under study.Precipitation during the recession period results in various shapes of the recession curves for the same spring.Therefore, analysis of a single recession period on a spring is not the most appropriate method for drawing major conclusions regarding the aquifer and groundwater supplies.
To overcome this problem the modified "matching strip" meth od was used (POSAVEC et al., 2006) for recession analy sis on the monitored springs.This method includes analysis of all recession events of the runoff hydrograph, allowing construction of a main recession curve.In some cases, curves with a lower coefficient of determination describe the recession more appropriately, and this can be easily determined visually (RIGGS, 1968).So the main recession curve, which has the highest coefficient of determination, is not always the best selection.Therefore, the final regression model is selected on the basis of both, experience and the existing software.
Recessions with two, three and sometimes even more discharge micro-regimes are quite common in karst aquifers (BONACCI, 1987).Such recession curves can be easily separated.For this purpose software was used for the automatic separation of the main recession curve of the monitored springs (POSAVEC et al., 2010;PARLOV, 2010).The programme is based on the analysis of the duration curves of mean daily discharge, i.e. on determining the critical discharge where the slope of the curve changes abruptly.
Analyses of complete hydrographs consist of statistical analysis of time series data (BOX & JENKINS, 1974).This method was introduced by MANGIN (1981MANGIN ( , 1984) ) to the karst hydrology and subsequently improved by PADILLA et al. (1994).Numerous authors have used this method for defining karst systems (LAROCQUE et al., 1997(LAROCQUE et al., , 2000;;TO-MASKO et al, 2001;PANAGOPOULOS & LAMBRAKIS, 2006;FIORILLO & DOGLIONI, 2010;STROJ, 2010;TER-ZIĆ et al., 2011)., Autocorrelation and crosscorrelation meth ods are used for the time series analysis of the springs, which is the aim of the current research (DAVIS, 2002;BOR-RADAILE, 2003).The characteristics of the springs catchment, the relationships between springs and the dependence between the precipitation in the catchment and the spring discharge were all determined using these methods.
The autocorrelation function represents the linear relationship of the successive data values inside a time series, dependant on their time distance (BOX & JENKINS, 1974;TERZIĆ et al., 2011).The autocorrelation method is based on the intercorrelation of data from a time series, with a successive increase of time distances, whereby the degree of similarity between the data for each step is calculated.Autocorrelation for the distance τ corresponds to the covariance of all measurements x t and the measurements with a time distance x t + τ according to the equation: (2) where cov τ is the covariance for a time distance, x is the time series, n is the number of measurements in a time series, τ is the time distance between two measurements, is the average value of the sample.
The autocorrelation coefficient (r x ) ranges from a maximum of 1 to the minimum of -1.The autocorrelation value of 1 means that the compared time series are identical.The time needed for the decrease of the correlation coefficient below 0.2 is called the "memory effect" MANGIN (1984), and reflects the duration of the reaction of the system to the input signal.Interpretation of the autocorrelogram should be undertaken with care, because the shape of the autocorrelogram depends, not only on the characteristics of the karst system, but also on the duration, frequency and the intensity of the rain (GRASSO & JEANINE, 1994;EISENLOHR et al., 1997).
To establish the places of pronounced similarities i.e. similarities between individual data, it is possible to compare two different time series in the same way one time series is compared with itself for a particular time lag.From such comparisons it is possible to obtain two kinds of information: the strength of the connection between these two time series and the time delay between them at the places of maximum similarity.The cross-correlation equation is similar to the autocorrelation coefficient equation.If two time series are marked as variables X and Y, and n is the number of pairs which are compared in one step (k) of the cross-correlation, than the crosscorrelation coefficient is as follows: (3)

METHODs UsED FOr DEFINING THE HYDrOGEOCHEMICAL CHArACTErIsTICs OF sPrING WATErs
Between February 2008 and November 2010, water samples were collected from four springs (Tonković, Pećina, Majerovo Vrilo and Klanac).Prior to sampling, the following parameters of the sample waters were measured ''in situ'' by probes from the WTW company: electrolytic conductivity (EC), temperature (T), pH and oxygen content.Water samples (500ml) were collected for measurements of basic cations, anions and alkalinity.Samples (100ml), were also collected for measurements of stable isotopes (δ 18 O and δD).The samples were collected under different hydrological con ditions, with a total of 16 samples per spring for basic cations and 25 samples for stable isotope measurements.Basic cations and anions were measured in the Hydrogeochemical laboratory, Department of Hydrogeology and Engineering Geology of the Croatian Geological Survey.Con centrations of sodium, potassium, magnesium and calcium were measured using the atomic absorber type AAS Analyst 700 by Perkin Elmer.Nitrate, orthophosphate, sulphate and chloride ions were measured using the ion chromatograph by LabAlliance.Bicarbonate ions were measured by titration.From the chemi cal composition of the water, interactions between the rock-mass and groundwater can be determined, as can the underground retention time of the water.Besides the natural processes, the water composition is also influenced by pollution, especially by nitrates found in artificial and natural fertilizers.Understanding of the reactions between the water and environmental substances introduced to it allows anticipation of the dangers/harmful influences of pollutants on water quality, and hence the influence on the human health and/or the health of the ecosystem.
The ratios of the stable isotopes deuterium (dD) and oxygen (d 18 O) in the sampled waters were measured using the mass spectrometar by Finnigan DELTA plus XP i DELTA plus in the JOANNEUM RESEARCH RESOURCES, Institute of Water, Energy and Sustainability, Department of Isotope Hydrology and Environmental Analytics in Graz (Austria) and in the Physics department of the School of Medicine in Rijeka (Croatia).These ratios are used in hydrogeological research to determine the recharge area and the hydrodynamic conditions, i.e. changes in water velocity in discrete parts of the aquifer and the hydraulic connection between separate aquifers (DOCTOR et al., 2000;VANDENSCHRICKA et al., 2002).In nature the differences in the ratio can be seen in the phase transitions in the water media.Therefore, variations of the isotope ratios in the groundwater are the result of seasonal changes, the altitude of the recharge area, proximity to the sea, intensive rain showers and isotope exchange with the rock-mass.The linear dependence between the stable isotopes oxygen-18 and deuterium in precipitation was recognized by FRIEDMAN (1953), while CRAIG (1961) defined this dependence with the following formula: which he named global meteoric water line (GMWL).The GMWL is calculated for the mean ratios of oxygen and hydrogen isotopes from precipitation around the world.However, the oxygen and hydrogen isotope ratios of some regions can differ from the global ones due to different meteorological conditions in the region.In this case, the local meteoric water line (LMWL) can be calculated.The physico-chemical indicators which were continuously collected during fieldwork, as well as hydrochemical measurements, were geochemically modelled using the NET PATH computer programme and statistical multivariate Factor Analysis using the STATISTICA software package.NETPATH was developed by PLUMMER et al. (1994) from the USGS and is based on the geochemical model of mass balance.It consists of a series of smaller programmes which, when put together, enable the input and management of chemical and isotope data from the analysed water.Comparison of the determinants between samples was undertaken by applica tion of the R-mode Factor Analysis, using the STATISTI CA 6.0 programme package (1998), and showed geochemical processes were the dominating ones determining the che mical composition of the water,.Factor analysis was performed on the following indicators: T, pH, Ca 2+ , Mg 2+ , nMg/nCa, Na+, K+, HCO 3 -, Cl -, SO 4 2-, NO 3 -, PO 4 3--P, SI cal , SI dol , logpCO 2.

recession analyses
The recession coefficients of the monitored springs were obtained using the "matching strip" method for the analysis of the falling limbs of the runoff hydrographs (Fig. 3).Since only the exponential regression model enables evaluation of the recession coefficient on which comparison of the properties of the aquifer are based,, this model was selected as being the most appropriate for comparison of the characteristics of the catchments of the monitored springs.For the discharge time series measured at Majerovo Vrilo, the most appropriate solution for the recession analysis was dividing the main recession curve into two parts at the point where the critical discharge rate was 6,4 m 3 /s.The value of the recession coefficient of the first part of the curve is α 1 = 0,247 day -1 , and the second α 2 = 0,029 day -1 .For other springs, division using three discharge micro-regimes was selected.At the Klanac spring the discharge regime changes at the value of 4,13 m 3 /s.This implies that high waters are best described by the recession coefficient of α 1 = 0,044.Changes in the runoff regime also appear at the critical discharge value of 1.14 m 3 /s.Hence the recession coefficient of α 2 = 0,032 was calculated for the medium waters.The recession coefficient of low waters is α 2 = 0,046.In the case of the Pećina spring, the main recession curve was divided at the critical discharge value of 6.00 m 3 /s, and the high waters are best described by the recession coefficient of α =0,115.The medium waters recession coefficient value is α 2 = 0,047 lasting until the critical discharge value of 0,65 m 3 /s.The value of the recession coefficient for low waters is α 3 = 0,07.Discharge of the aquifer at the Pećina spring during the middle water regime is not uniform (Fig. 3), suggesting that the regression model selected on the coefficient of determination is not significant.At the Tonković spring, changes in the discharge regime appear at values of 4.20 m 3 /s and 2,70 m 3 /s.The recession coefficient defined by the automatic recession analysis is α 1 = 0,023 for high waters, α 2 = 0,014 for medium waters and α 3 = 0,009 for low waters.As at the Pećina spring, during the high water regime, significant differences can be determined in the drainage speed during individual recession periods at the Tonković spring.This was the reason why it was not possible to analyse the high water recession curve with more certainty.Recession analysis (Fig. 3, Table 1) clearly showed that the Pećina and Majerovo Vrilo springs have significant differences between the recession coefficients of the quick and the base flow, when compared to the Klanac and Tonković springs.The highest runoff differences between the base and quick flow occurs at the Majerovo Vrilo spring with the difference in the recession coefficient between the base and the quick flow being of one order of magnitude.Furthermore, as seen on the recession curve graphs of the springs, the Peći na and Tonković springs have significant differences in the types of drainage of the fissured systems during different quick flow recession events.At the Pećina spring, during individual recession periods, the runoff decreased from 5 m 3 /s to 1 m 3 /s in only 4 -5 days, while during other recession periods for the same runoff more days are needed.At the Tonković spring, significant differences in the discharge rates during different recession periods of the high water regimes can also be determined.While in some cases of recession, runoff decreases from 7 to 4 m 3 /s in only 2 days, sometimes 20 days were needed for the same decrease.This variation in the runoff properties at the same spring could be the result of: the varying intensity and spatial distribution of precipitation in the catchment, the land cover, differences in evapo transpiration rates, and the heterogeneity of the aquifer, which directly influences the dynamics of the runoff.At the Klanac and Majerovo Vrilo springs, no such significant differences in the runoff during different recession events were determined.It was established, after the statistical analysis of the discharge time series using the autocorrelation method, that all four investigated springs drain water from aquifers with similar hydrogeological properties.This is evident from the similarities of their autocorrelation functions (AFK) (Fig. 4).Nevertheless, at the Pećina and Majerovo Vrilo springs the AKF have a steeper slope at the beginning and the values of the memory effect are lower than at the Klanac and Tonković springs.This is the result of the better developed underground karst channel network of the Pećina and Majerovo Vrilo aquifers, which is in accordance with the results of the recession analysis.When interpreting AKF, it should be known that they relate to the average properties of the analyzed time series.
From the cross-correlation functions (KKF) of discharge time series of the different springs, the connection and the interrelation of the springs can be recognized (Fig. 5).Symmetrical KKF values of all springs indicate that the springs are independent of each other and depend on some other variable.In this case it refers mostly to precipitation in the catch ment of the monitored springs, because their quantity and time and space distribution, together with the hydrogeological aquifer properties, influence the runoff dynamics at the springs.It is supposed, based on the high values of the maximum crosscorrelation coefficients (max r xy ), which range from 0,77 -0,91, that all springs drain aquifers from catchments of similar hydrogeological characteristics and In cases of extreme differences in the recession characteristics of the aquifer, it is not possible to determine the recession curve with a high value of the coefficient of determination R 2 .Nevertheless, the "matching strip" method is still considered to be the more appropriate solution for the determination of the recession properties of the aquifer compared to analysing individual recession events.For the analysis of individual recession events, it is not possible to determine how much the runoff regimes of different recession events at the same spring differ from each other.

Analysis of the complete hydrographs
From the basic statistical water flow parameters presented in Table 2 it can be seen that the mean discharge values (Q) of the Klanac, Tonković and Majerovo Vrilo springs are very si milar, while at the Pećina spring it is almost twice as low.The values of the variation coefficient (c v ), which represents the ratio between the standard deviation (σ) and the mean value ( ) of the sample, indicate that the highest discharge variations occur at the Klanac, and lowest at the Tonković springs.
with a similar precipitation volume.The delayed response to rainfall of the Klanac and Tonković springs if compared to the Majerovo Vrilo spring, ranges from 24 -36 hours.The delay at the Pećina spring is between 12 -24 hours.The Majerovo Vrilo spring reacts up to 12 hours prior to the Pećina spring, the Tonković spring up to 12 hours prior to the Klanac spring.The highest crosscorrelation coefficients are those of the Klanac and Tonković springs, with the lowest at Pećina and Majerovo Vrilo, which is expected since these two springs have the most distant catchments.
The maximum discharge and precipitation cross-correlation coefficients max r xy (Fig. 6) are not high (0.3 -0.46).The reason for the low values is that the discharge dynamics at the springs does not solely depend on the quantity of precipitation in the catchment.It also depends on the season, i.e. air temperature during precipitation, the vegetation and the hydrogeological properties, including the degree of karstification, the existence of cover deposits, their thickness and other factors.Nevertheless, despite the low cross-correlation coefficients the connection between discharge and precipitation is clearly pronounced.At the Pećina spring, the maximum crosscorrelation coefficient was registered during the three day delay.This means that the average reaction delay of this spring is three days.At the Majerovo Vrilo spring, the reaction to precipitation is the same as for the Pećina spring, but the max r xy is somewhat higher than the other springs.The crosscorrelation coefficients of the Klanac and Tonković springs have two maxima which means that different precipitation events cause variation in delay.Hence the Klanac spring reacts in 4 to 10 days, and Tonković spring 5 to 11 days after the onset of precipitation.Various delay times at these springs could be the result of variation in infiltration of precipitation in the karst underground over both time and place in the catchment.Rapid reaction of a spring to precipitation indicates fast infiltration in the epikarst zone and further into the underground channel system, which enables rapid transport.Another reason could be a smaller catchment surface which makes the path the infiltrated water must travel from the surface to the spring significantly shorter.After a time period of 40 days the function gradually starts to grow again as seen on the cross-correlograms.This is in accordance with the autocorrelation functions of the springs, i.e. it is the result of repeated incoming water waves, especially in the wetter part of the year.
All hydrograph statistical analysis, as well as recession hydrograph analysis used in this study, indicate faster reactions to rainfall by the Majerovo Vrilo and Pećina compared to the Tonković and Klanac springs.This implies a better developed underground channel network of the aquifer systems of the Majerovo Vrilo and Pećina springs.The Klanac  and Tonković springs aquifer systems show better retention properties and lower permeability due to the network of smal ler fractures within the aquifer.

The isotopic composition of the waters
The LMWL was calculated for the research area by MANDIĆ et al. (2008) (Fig. 7) who emphasised that there is a need for several more years of measurements in order to determine a more precise value.The LMWL equation is as follows: dD = 6.8·d 18 O + 1.5 (5) The measured stable isotope oxygen-18 and deuterium values of the spring waters are scattered around the LMWL and GMWL values, clearly indicating that in the study area the groundwater is recharged from precipitation (Fig. 7).It was shown by MANDIĆ et al. (2008) that the average deuterium excess in the investigated area is about 12 ‰.Such values are evidence of the influence of the Mediterranean type of climate, especially during the winter months, but the main influence is from the continental type of climate.
No significant fluctuations in oxygen18 values (Fig. 8) were observed in the spring waters of the river Gacka catchment which is in agreement with the observations of MAN-DIĆ et al. (2008).The lowest negative d 18 O values were mea sured at Majerovo Vrilo, followed by those at the Tonko vić spring.The Pećina spring has the highest positive values.These values imply the highest altitude of the recharge area occurs at the Majerovo Vrilo and the lowest at the Pećina spring.For Europe and Africa, the δ 18 O isotope gradient ranges from -0.15 to -0.5 ‰ for every 100 m increase in altitude (CLARK & FRITZ, 1997).In order to calculate the altitude of the recharge area, the isotope gradient of -0,25 ‰ was selected and the δ 18 O values measured in spring and autumn, as suggested by CLARK & FRITZ (1997).The altitude differences obtained range from 140 to 260 m.CLARK & FRITZ, (1997), determined that the depletion for 18 O varies between about -0.15 to -0.5 ‰ per 100 m rise in altitude.The calculated altitude values of the dominant recharge area of the springs are as follows: • Majerovo Vrilo -1020 m a.s.l., • Tonković -760 m a.s.l.• Pećina -620 m a.s.l.

The quality of the spring waters
The results of the measurements of the water quality indicators are given in Table 3.
Based on their basic chemical composition the spring waters of Majerovo Vrilo belong to the CaHCO 3 to CaMg-HCO 3 -type of water, i.e. to the hydrochemical facies (Fig. 9).Spring waters of Tonković, Klanac and Pećina belong to the Ca-HCO 3 type of water.These water types are the result of the dissolution of carbonate minerals (calcite and dolomite) which are the main constituents of the catchment area of the springs.The temperature of the spring waters, was measured periodically during the research period, and varies from 8.4 to 10.0 °C.This is generally in accordance with the mean annual air temperature of the recharge area of the spring.The lowest average measured temperatures were recorded at Majerovo Vrilo, which is in accordance with the data suggesting that this spring is recharged from the highest altitudes, where the temperatures are the lowest.Temperatures measured at the Tonković and Pećina springs are also compatible with the calculated altitudes of their recharge area.
The oxygen content measured in the spring waters, as well as the model determined carbonate saturation level, is advantageous, from the water quality protection standing point.A high oxygen concentration was measured in the waters which enables oxidation processes useful as water autopurification reactions.Spring waters are mostly saturated with calcite, which enables deposition of any existing metals in the water due to their affinity with carbonates.Hence, in saturated conditions these reactions could be initiated, which also has a water autopurification role.
The measured nitrate concentrations are below the Croatian drinking water standard.(OFFICIAL GAZETTE NO. 47/08) and range from 1.3 to 5.2 mg/l with small seasonal oscillations.The highest concentration was measured at the Pećina spring (Table 3) in the autumn and winter.This is the consequence of both anthropogenic influences and degradation of the plant material in the highly karstified aquifer of this spring.
Ammonium and orthophosphate concentrations are very low in all the springs and sometimes they are below the detection limit.Sulphate concentrations range from 3.4 -12.9 mg/l and are below levels for the Croatian drinking water standard (OFFICIAL GAZETTE NO. 47/08).Concentrations of chloride range from 1.3 -14.1 mg/l with the highest concentration being observed at the Tonković spring in March (14.1 mg/l) and April (8.9 mg/l).Furthermore, concentrations of chloride at other springs are higher in March and April, than during other months.These high concentrations are most like ly to be a consequence of de-icing of the roads, because the springs are located near them.
Waters of the river Gacka springs are still characterised by their excellent quality and are classified as strategically important water reserves of the Republic of Croatia (OFFI-CIAL GAZETTE NO. 91/08).Hence, it is important to monitor and to promptly react to existing and potential pollution sources in the catchment.Taking into account previous hydrogeochemical studies carried out in this area (LUKAČ RE BERSKI, 2008;LUKAČ REBERSKI et al., 2009), it is evident that, despite the good protection and quality, the analysed water quality indicators show an upward trend due to the anthropogenic influence in the spring catchments.

statistical analysis of the geochemical indicators of spring waters
Five factors were selected based on the results of the Factor Analysis i.e. the factor loadings (Table 4).The first is the lithogenic factor and refers to the influence of the partial pres sure on the calcite and dolomite saturation indices, i.e. the influence on the dissolution/precipitation of limestone and dolomite.The second factor is also a lithogenic factor, and represents the influence of dolomite dissolution on the water chemistry.The third factor is anthropogenic; it shows the influence of the wash out of the terrain and epikarst zone on the chemistry of the water.The fourth factor is lithogenic and represents the dissolution of limestone.The fifth factor is anthropogenic and indicates the influence of salting of roads on the chemical composition of the water.

The calculations of mixing ratios
Mixing models for individual springs were provided using the NETPATH programme.The water ratios from the Majerovo Vrilo and Tonković sprin catchments in the Klanac spring were determined in order to establish if this spring shares the catchment area with the others.The water ratios vary depending on the hydrological conditions.During high and medium waters the Klanac spring mostly shares the catchment area with the Tonković spring (Fig. 10).During low water periods it shares its catchment area equally with Majerovo Vrilo and the Tonković springs.It should be noted that these values are valid for the time of sampling.For a more detailed analysis of their interrelationships samples should be collected on a daily basis throughout the hydrological year.

The subdivision of the catchment
The catchment of the river Gacka springs was divided into three subcatchments based on the results of the investigations, the interpretation of the lithological, structural and hydrogeological characteristics of the aquifer system and the hydrogeochemical properties of the spring waters, (Fig. 11).
The Tonković and Klanac springs, which are about 50 m away from each other, drain the same catchment as shown by the research results.Hence, it was not possible to establish the groundwater divide between these two springs.The   LOGOVIĆ, 1989), it is difficult to determine the cause of this discrepan with a high degree of certainty.There are several possible explanations.One could be that during periods of heavy precipitation, the amount of rainfall exceeds the infiltration rate and the water level rises in the epikarst  A common zone was also defined between the Majerovo Vrilo catchment and the Klanac and Tonković springs.Water tracings performed in Vrhovinsko polje and at Trnavac confirmed an underground connection common to both catch ments.Furthermore, the hydrogeochemical mixing model showed that during high water levels the groundwater flows to the Klanac spring from the same part of the catchment as to the Tonković spring.Under low water conditions, at the Klanac spring there is an increased ratio of groundwater from the Majerovo Vrilo subcatchment.This is probably due to hydraulic conditions which prevent water recharge from the Majerovo Vrilo catchment during the high water regime.When the water level drops, water recharge from Majerovo Vrilo (Fig. 12) is enabled.It is supposed that the water from the Klanac catchment does not discharge towards the Majerovo Vrilo.

CONCLUsION
Applied research techniques have enabled the recognition of subcatchments in, parts of the large karst aquifer system of the river Gacka springs and the determination of their dynamic characteristics, as well as their interference.All the former hydrological analyses were performed based on hydrological da measured only on the river Gacka profiles.Such data were referred to the entire system, while the characteristics of certain parts of the catchment still remain unknown.Therefore, in the cases of such large catchment areas, it is very useful to divide it into smaller units or subcatchments.Results presented here show that there are significant differences in hydrogeological characteristics of individual parts of the large river Gacka springs catchment area.Based on interpretation of the lithological and structural-tectonic characteristics of the area, the tracing test data, and the applied research methods, the river Gacka catchment is divided into three subcatchments: the Majerovo Vrilo, the Tonković and Klanac, and the Pećine subcatchments.The boundaries between them are not presented as lines, but were determin ed as zones which are common to the neighbouring subcatchments (Fig. 11).
Statistical analysis, analyses of the recession parts of the hydrograph, together with hydrogeochemical analyses, all indicate a better development of the groundwater flow paths of the aquifer systems in the Majerovo Vrilo and Pećina springs catchments when compared to the Tonković and Klanac springs.This is the reason why the Majerovo Vrilo and Pećina spring catchments are more sensitive to pollution.Due to their well developed network of underground water flow and the short subsurface residence time of the water, they have a low water autopurification capacity.The Tonković spring during the quick flows, showed significant variation in drainage during different recession events, as observed from recession analysis.This led to the conclusion that parts of this subcatchment are even more sensitive to pollution.
The autocorrelati and cross-correlation time series analysis and the matching strip method are very useful methods for the hydrogeological research of karst aquifers.They enable collection of data regarding the average aquifer characteristics which can be used f interpretation of the hydrogeological properties of the catchment, such as the structural properties, type of porosity, development of the channel network, dynamics of the aquifer system, sensitivity to pollu- tion and others.Analysis of the recession hydrograph using the matching strip method proved to be very useful because, in addition to calculation of the recession coefficient, the diagram can be used to establish the time period of various recession events, since all recession events of the time series are presented on the diagram.Hence, such recession analysis is more reliable than cases where several recession events are considered individually.
The water of the river Gacka springs is still of excellent quality.Hydrogeochemical investigations imply a trend of growing concentrations of anthropogenic indicators (nitrates, orthophosphates, sulphates and chlorides) in the spring's catch ments.Furthermore, the dominant processes which influence the composition of the spring waters were determin ed, together with the average recharge altitudes of the springs, and using a geochemical mixing model, the interference between the catchments was also established.
Results presented here point out the necessity for a multidisciplinary research approach, especially in order to determine the hydrogeological characteristics of complex aquifer systems like the river Gacka springs catchment.Due to the high importance of this catchment, which was proclaimed to be an area of strategic drinking water reserves for the Republic of Croatia, the collected data are useful for all future development plans, both for groundwater exploitation as well as for protection from possible contamination.

Figure 3 :
Figure 3: Recession analyses at the observed springs.
At the Tonković spring the autocorrelation function (AKF) has a steeper slope for the first 20 days.The memory effect, i.e. the average duration of the springs reaction to precipitation is 62 days.Similarly, the slope of the function is steeper during the first 22 days at the Klanac spring and the memory effect is 63 days.At the Majerovo Vrilo and Pećina springs, the shape of the autocorrelation functions are somewhat different.At the Majerovo Vrilo spring the steepest slo pe of the function is in the first 4 days, while the response of the catchment, from which the spring is recharged, is 26 days.At the Pećina spring the fast discharge lasts 9 days on average, and the memory effect of the spring is 33 days.Steep slopes of the autocorrelation function (AKF) indicate rapid infiltration of precipitation and faster drainage of the underground system through the well developed network of karst channels.Its more gentle slopes reflect slower drainage of the fissured porosity aquifer area where discharge is similar to the intergranular aquifers.

Figure 6 :
Figure 6: Crosscorellograms of springs which represent the relationship between rainfall and discharge.

Figure 7 :
Figure 7:The ratio of oxygen-18 and deuterium stable isotopes in spring waters of the Gacka catchment.

Figure 8 :
Figure 8: The distribution of stable isotope oxygen-18 in the spring waters.

Figure 10 :
Figure 10: Plot of the calculated mixing ratios of waters from the Tonković and Majerovo Vrilo springs in the water from Klanac spring.

Figure 12 :
Figure 12: A conceptual model of the river Gacka springs catchment area.

Table 1 :
The values of the recession coefficients at observed springs.

Table 2 :
The values of basic statistical parameters for the flow rate.

Table 3 :
Results of measurements of hydrochemical indicators.
groundwater divide between the catchment of the Pećina spring and the Tonković and Klanac springs is not a line, but represents a recharge zone which is determined between the se springs and is common to both catchments.Most of the Tonković and Pećina springs catchment area drains the Tonković spring.During high water levels a connection exists between the sink in Kozjan with the Pećina spring, as evidenced by tracing test results.Furthermore, these test results indicate that the Pećina spring has a larger catchment surface than expected based on the mean annual discharge values.Considering the complex structural-tectonic relationships, which are a reflection of many phases of karstification and tectonism, particularly the neo tectonic dynamics (PRE-