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

Headnote

ABSTRACT

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.

Keywords: hydrogeology, karst spring, hydrology analysis, hydrogeochemistry, river Gacka springs

(ProQuest: ... denotes formulae omitted.)

1. INTRODUCTION

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 (Tonkovic, Majerovo Vrilo, Klanac and Pecina) 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^sup 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 (PAVICIC et al., 2003; LUKAC REBERSKI, 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 & Andric, 2008; Lukac 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 (Tonkovic, Majerovo Vrilo, Klanac and Pecina), 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 (LUKAC REBERSKI, 2011).

2. 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 Velebit, Kapela, and Pljesivica 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, (ZANINOVIC et al., 2004) while the annual mean precipitation values range between 1000 and 2500 mm (GAJIC-CAPKA 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).

3. 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 (VLAHOVIC 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 catchment. 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 (LUKAC REBERSKI 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 position changes depending on the groundwater level (ZUGAJ, 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 Na-fluorescein 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 Covici profile extends to 516 km^sup 2^ (BIONDIC et al. 2010).

Groundwater-flow 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 (PRELOGOVIC, 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).

4. 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, whereas 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 m3s-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 (LUKAC 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 annual precipitation was calculated using the regression equation defined for the research area (GAJIC-CAPKA 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 Covici-Podgora is 14.2 m3/s.

The discharge coefficient c can be defined as follows:

... (1)

where Q is the mean discharge (m3s-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^sup 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.

5. 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 extremely heterogeneous characteristics of the karst aquifer, spring runoffhydrographs 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 runoffhydrographs. 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 (TALLAKSEN, 1995; DEWANDEL et al., 2003; BAKALOWICZ, 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 (KRESIC, 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" method was used (POSAVEC et al., 2006) for recession analysis on the monitored springs. This method includes analysis of all recession events of the runoffhydrograph, 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 (1981, 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, 2000; TOMASKO et al, 2001; PANAGOPOULOS & LAMBRAKIS, 2006; FIORILLO & DOGLIONI, 2010; STROJ, 2010; TERZIC et al., 2011)., Autocorrelation and cross-correlation methods are used for the time series analysis of the springs, which is the aim of the current research (DAVIS, 2002; BORRADAILE, 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; TERZIC 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 xt and the measurements with a time distance xt + τ 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 (rx) 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 cross-correlation coefficient is as follows:

... (3)

6. METHODS USED FOR DEFINING THE HYDROGEOCHEMICAL CHARACTERISTICS OF SPRING WATERS

Between February 2008 and November 2010, water samples were collected from four springs (Tonkovic, Pecina, 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 (δ18O and δD). The samples were collected under different hydrological conditions, 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. Concentrations 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 chemical 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 (d18O) 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:

... (4)

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 NETPATH 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 application of the R-mode Factor Analysis, using the STATISTICA 6.0 programme package (1998), and showed geochemical processes were the dominating ones determining the chemical composition of the water,. Factor analysis was performed on the following indicators: T, pH, Ca2+, Mg2+, nMg/nCa, Na+, K+, HCO3 -, Cl-, SO4 2-, NO3 -, PO4 3-P, SIcal, SIdol, logpCO2.

7. RESULTS AND DISCUSSION

7.1. 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 runoffhydrographs (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 m3/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 m3/s. This implies that high waters are best described by the recession coefficient of α1 = 0,044. Changes in the runoffregime also appear at the critical discharge value of 1.14 m3/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 Pecina spring, the main recession curve was divided at the critical discharge value of 6.00 m3/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 m3/s. The value of the recession coefficient for low waters is α3 = 0,07. Discharge of the aquifer at the Pecina 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 Tonkovic spring, changes in the discharge regime appear at values of 4.20 m3/s and 2,70 m3/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 Pecina spring, during the high water regime, significant differences can be determined in the drainage speed during individual recession periods at the Tonkovic 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 Pecina and Majerovo Vrilo springs have significant differences between the recession coefficients of the quick and the base flow, when compared to the Klanac and Tonkovic springs. The highest runoffdifferences 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 Pecina and Tonkovic springs have significant differences in the types of drainage of the fissured systems during different quick flow recession events. At the Pecina spring, during individual recession periods, the runoffdecreased from 5 m3/s to 1 m3/s in only 4 - 5 days, while during other recession periods for the same runoffmore days are needed. At the Tonkovic 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, runoffdecreases from 7 to 4 m3/s in only 2 days, sometimes 20 days were needed for the same decrease. This variation in the runoffproperties 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 evapotranspiration 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 runoffduring different recession events were determined.

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 R2. 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 runoffregimes of different recession events at the same spring differ from each other.

7.2 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 similar, while at the Pe.ina spring it is almost twice as low. The values of the variation coefficient (cv), 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.

At the Tonkovic 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 Pecina springs, the shape of the autocorrelation functions are somewhat different. At the Majerovo Vrilo spring the steepest slope 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 Pecina 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.

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 Pecina 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 Tonkovic springs. This is the result of the better developed underground karst channel network of the Pecina 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 catchment of the monitored springs, because their quantity and time and space distribution, together with the hydrogeological aquifer properties, influence the runoffdynamics at the springs. It is supposed, based on the high values of the maximum cross-correlation coefficients (max rxy), which range from 0,77 - 0,91, that all springs drain aquifers from catchments of similar hydrogeological characteristics and with a similar precipitation volume. The delayed response to rainfall of the Klanac and Tonkovic springs if compared to the Majerovo Vrilo spring, ranges from 24 - 36 hours. The delay at the Pecina spring is between 12 - 24 hours. The Majerovo Vrilo spring reacts up to 12 hours prior to the Pecina spring, the Tonkovic spring up to 12 hours prior to the Klanac spring. The highest cross-correlation coefficients are those of the Klanac and Tonkovic springs, with the lowest at Pecina and Majerovo Vrilo, which is expected since these two springs have the most distant catchments.

The maximum discharge and precipitation cross-correlation coefficients max rxy (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 Pecina spring, the maximum cross-correlation 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 Pecina spring, but the max rxy is somewhat higher than the other springs. The cross-correlation coefficients of the Klanac and Tonkovic 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 Tonkovic 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 Pecina compared to the Tonkovic and Klanac springs. This implies a better developed underground channel network of the aquifer systems of the Majerovo Vrilo and Pecina springs. The Klanac and Tonkovi. springs aquifer systems show better retention properties and lower permeability due to the network of smaller fractures within the aquifer.

7.3 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:

... (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 oxygen-18 values (Fig. 8) were observed in the spring waters of the river Gacka catchment which is in agreement with the observations of MANDI. et al. (2008). The lowest negative d18O values were measured at Majerovo Vrilo, followed by those at the Tonkovi. 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 Â18O 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 Â18O 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 18O 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.

7.4. 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 Ca-HCO3 to CaMg- HCO3-type of water, i.e. to the hydrochemical facies (Fig. 9). Spring waters of Tonkovi., Klanac and Pe.ina belong to the Ca-HCO3 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 Tonkovic 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 likely 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 (OFFICIAL 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 (LUKAC REBERSKI, 2008; LUKAC 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.

7.5. 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 pressure 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.

7.6. The calculations of mixing ratios

Mixing models for individual springs were provided using the NETPATH programme. The water ratios from the Majerovo Vrilo and Tonkovic 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 Tonkovic spring (Fig. 10). During low water periods it shares its catchment area equally with Majerovo Vrilo and the Tonkovic 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.

7.7. 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 Tonkovic 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 groundwater divide between the catchment of the Pecina spring and the Tonkovic and Klanac springs is not a line, but represents a recharge zone which is determined between these springs and is common to both catchments. Most of the Tonkovic and Pecina springs catchment area drains the Tonkovic spring. During high water levels a connection exists between the sink in Kozjan with the Pecina spring, as evidenced by tracing test results. Furthermore, these test results indicate that the Pecina 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 (PRELOGOVIC, 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 zone. This causes water overflow into the preferred flows near the surface, and possibly establishing a connection with the Pecina spring, probably due to the infiltration through the shallower parts of the aquifer. Another reason could be water flowing over a supposed hydrogeological barrier in the spring area during high waters. Furthermore, underground discharge at the Pecina spring is not excluded due to the fractured nature of the area surrounding the spring. This could be the case especially during low waters and could occur without any obvious surface runoff. In that case a larger catchment surface is drained without any visible evidence on the surface. This could also be in accordance with the recently postulated theory (BONACCI & ANDRIC, 2008) which assumed that part of the waters of the river Lika disappears and flows underground to the river Gacka catchment, draining downstream in the Pecina spring. Three permanent springs are located downstream of the Tonkovic spring: Jaz, Marusino Vrilo and Graba. They are characterised by small fluctuations in discharge under different hydrological conditions which indicate the existence of a barrier in their hinterland which slows down the groundwater and prevents significant water level fluctuations.

A common zone was also defined between the Majerovo Vrilo catchment and the Klanac and Tonkovic springs. Water tracings performed in Vrhovinsko polje and at Trnavac confirmed an underground connection common to both catchments. 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 Tonkovic 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.

8. 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 Tonkovic and Klanac, and the Pecine subcatchments. The boundaries between them are not presented as lines, but were determined 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 Pecina springs catchments when compared to the Tonkovic and Klanac springs. This is the reason why the Majerovo Vrilo and Pecina 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 Tonkovic 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 catchments. Furthermore, the dominant processes which influence the composition of the spring waters were determined, 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.

References

REFERENCES

BAHUN, S. (1974): Tektogeneza Velebita i postanak Jelar naslaga [The tectogenesis of the Mt. Velebit and the formation of Jelar deposits - in Croatian].- Geol. vjesnik, 27, 35-51.

BAHUN, S. & FRITZ, F. (1975): Hidrogeoloske specificnosti Jelarnaslaga [Hydrogeological specificities of the Jelar breccias - in Croatian].- Geol. vjesnik, 28, 345-355.

BAKALOWICZ, M. (2005): Karst groundwater: a challenge for new resources.- Hydrogeology Journal 13, 148-160. doi:10.1007/ s10040- 004-0402-9.

BIONDIC, R., BIONDIC, B. & MEASKI, H. (2010): Novelacija granica zona sanitarne zastite izvorista Gacke [Gacka spring sanitary protection zone novelation - in Croatian]. Geotehnical faculty, University of Zagreb, 115 p.

BONACCI, O. (1987): Karst Hydrology with Special Reference to the Dinaric Karst.- Springer-Verlag-Berlin-Heidelberg, 9/4, 328-338.

BONACCI, O. (1995): Meðuzavisnost geologije, hidrogeologije i hidrologije pri rjesavanju hidrotehnickih problema. [The interdependence of geology, hydrogeology and hydrology in solving hydraulic proble - in Croatian].- First Croatian geological congress, Book of Proceedings, 1, 105-108.

BONACCI, O. & ANDRIC, I. (2008): Sinking karst rivers hydrology: case of the Lika and Gacka/Croatia.- Acta Carsologica, 37/2, 185- 196.

BOX, G.E.P. & JENKINS, G.M. (1974): Time Series Analysis: Forecasting and Control.-Holden Day, San Francisco, 575 p.

CLARK, I. & FRITZ, P. (1997): Environmental Isotopes in Hydrogeology.- CRC Press, New York, 20, 174-188.

CRAIG, H. (1961): Isotope variations in meteoric waters.- Science, 133, 1702-1703.

DAVIS, J.C. (2002): Statistics and Data Analysis in Geology. Cansas geological survey.- The University of Cansas, third edition, USA, 638 p.

DEWANDEL, B., LACHASSAGNEB, P., BAKALOWICZ, M., WENGB, P. & AL-MALKI, A. (2003): Evaluation of aquifer thickness by analysing recession hydrographs. Application to the Evaluation Oman ophiolite hard-rock aquifer.- J. Hydrol., 274, 248-269. doi: 10.1016/S0022-1694(02)00418-3.

DOCTOR, D. H., LOJEN, S. & HORVAT, M. (2000): A stable isotope investigation of the classical Karst aquifer: evaluating karst groundwater components for water quality preservation.- Acta carsologica, 29/1, 5, 79-92.

EISENLOHR, L., BOUZELBOUDJEN, M., KIRALY, L. & ROSSIER, Y. (1997): Numerical versus statistical modelling of natural response of a karst hydrogeological system.- J. Hydrol. 202, 244-262.

FIORILLO, F. & DOGLIONI, A. (2010): The relation between karst spring discharge and rainfall by cross-correlation analysis (Campania, southern Italy).- Hydrogeology Journal, 18/8, 1881-1895. doi:10.1007/ s10040-010-0666-1.

FRIEDMAN, I. (1953): Deuterium content of natural waters and other substances.- Geochim. et. Cosmochim. Acta, 4, 89-103.

GAJIC-CAPKA, M., PERCEC TADIC, M. & PATARCIC, M. (2003): Digitalna godisnja oborinska karta Hrvatske. [A Digital Annual Precipitation Map of Croatia - in Croatian] Hrvatski meteoroloski casopis (Croatian Meteorological Journal), 38, 21-33.

GRASSO, D.A. & JEANNIN, P-Y. (1994): Etude critique des méthodes d'analyse de la réponse globale des systèmes karstiques. Application au site de Bure (JU, Suisse).- Bulletin d'Hydrogéologie de l'Université de Neuchâtel, 13, 87-113.

KENDALL, C. & MCDONNELL, J. J. (1998): Isotope Tracers in Catchment Hydrology.-Elsevier Science B.V., Amsterdam, 839 p.

KRESIC, N. (1997): Quantitative solutions in hydrogeology and groundwater modeling.-Lewis Publishers, New York, 461 p.

LAROCQUE, M., MANGIN, A., RAZACK, M. & BANTON, O. (1997): Contribution of correlation and spectral analyses to the regional study of a large karst aquifer (Charente, France).- J. Hydrol., 205, 217-231. doi: 10.1016/S0022-1694(97)00155-8.

LAROCQUE, M., BANTON, O. & RAZACK, M. (2000): Transient State History Matching of a Karst Aquifer Ground Water Flow Model.- Ground Water, 38/6, 939-946.

LUKAC REBERSKI, J. (2008): Hidrogeoloska i hidrogeokemijska osnova za definiranje slijeva Gacke i zastita njenog izvorista [Hydro geological and hydrogeochemical basis for the catchment area definition: River Gacka springs catchment area - in Croatian].- Unpubl. Master's Thesis. Rudarsko-geolosko-naftni fakultet, Sveuciliste u Zagrebu (Faculty of mining, geology and petroleum engineering, University of Zagreb), 118 p.

LUKAC REBERSKI, J., KAPELJ, S. & TERZIC, J. (2009): An estimation of groundwater type and origin of the complex karst catchment using hydrological and hydrogeochemical parameters: A case study of the river Gacka springs.- Geol. Croat., 62/3, 157-178, doi: 10.4154/ gc.2009.15.

LUKAC REBERSKI, J. (2011): Odreðivanje podsljevova izvorista rijeke Gacke na osnovi hidrogeoloskih parametara [Definition of the river Gacka springs subcatchment areas on the basis of hydrogeological parametars - in Croatian].- Unpubl. Doctoral Thesis. Rudarskogeolosko- naftni fakultet, Sveuciliste u Zagrebu (Faculty of mining, geology and petroleum engineering, University of Zagreb), 158 p.

MANDIC, M., BOJIC, D., ROLLER-LUTZ, Z., LUTZ, H. & BRONIC KRAJCAR, I. (2008): Note on the spring region of Gacka River (Croatia).- Isotop. Envir. Health Studies, 44/2, 201-208. doi:10.1080/ 10256010802066364.

MANGIN, A. (1981): Utilisation des analyses corré1atoire et spectrale dans l'approche des systémes hydrologiques.- Computers Rendues Acad. Sci. Paris, 293, 401-404.

MANGIN, A. (1984): Pour une meilleure connaissance des systemes hydrologiques a partir des analyses correlatoire et spectrale.- J. Hydrol. 67, 25-43.

OFFICIAL GAZETTE NO. 91 (2008): Strategija upravljanja vodama [Water Management Strategy - in Croatian].

OFFICIAL GAZETTE NO. 47 (2008): Pravilnik o zdravstvenoj ispravnosti vode za pice [Drinking water regulations - in Croatian].

PADILLA, A., PULIDO-BOSCH, A. & MANGIN, A. (1994 Relative importance of Basefland Quickflow from Hydrographs of Karst Spring.- Ground Water, 32/2, 267-277.

PANAGOPOULOS, G. & LAMBRAKIS, N. (2006): The contribution of time series analysis to the study of the hydrodynamic characteristics of the karst systems: Application on two typical karst aquifers of Greece (Trifilia, Almyros Crete).- J. Hydrol., 329/3-4, 368-376. doi: 10.1016/j.jhydrol.2006.02.023.

PARLOV, J. (2010): Identifikacija parametara za modeliranje toka podzemne vode glavnih izvora u porjecju Mirne [Identification of Parameters Used for Groundwater Modeling of Major Springs in the Mirna River Valley - in Croatian].- Unpubl. Doctoral Thesis. Rudarsko- geolosko-naftni fakultet, Sveuciliste u Zagrebu (Faculty of mining, geology and petroleum engineering, University of Zagreb), 111 p.

PAVICIC, A., KAPELJ, S. & LUKAC, J. (2003): The influence of the Highway on the protected spring of Gacka river.- RMZ - Materials and geoenvironment, 50/1, 289-292.

PLUMMER, L.N., PRESTEMON, E.C. & PARKHURST, D.L. (1994): An interactive code (NETPATH) for modelling net geochemical reactions alog a flow path, Version 2.0, USGS.- Water-Resources Investigations Report 94-4169, Reston, Virginia.

POSAVEC, K., BACANI, A. & NAKIC, Z. (2006): A Visual Basic Spreadsheet Macro for Recession Curve Analysis.- Ground Water, 44/5, 764-767. doi: 10.1111/j.1745-6584.2006.00226.x.

POSAVEC, K., PARLOV, J. & NAKIC, Z. (2010): Fully Automated Objective- Based Method for Master Recession Curve Separation.- Ground Water, 48/4, 598-603. doi: 10.1111/j.1745-6584.2009. 00669.x.

PRELOGOVIC, E. (1989): Neotektonski pokreti u podrucju sjevernog Velebita i dijela Like [Neotectonic movements in northern Velebit and part of Lika - in Croatian]. Geol. vjesnik, 42, 133-147.

RIGGS, H.C. (1968): Techniques of Water-Resources Investigations of the United States Geological Survey, Chapter A1, Some statistical tools in hydrology. Book 4, 39 p.

STROJ, A. (2010): Podzemni tokovi u zaleðu krskih priobalnih izvora na podrucju Velebitskog kanala [Underground water flows in the hinterland of the velebit channel coastal karst springs - in Croatian]. Unpubl. Doctoral Thesis. Rudarsko-geolosko-naftni fakultet, Sveuciliste u Zagrebu (Faculty of mining, geology and petroleum engineering, University of Zagreb), 259 p.

TALLAKSEN, L.M. (1995): A review of baseflow recession analysis.- J. Hydrol., 165/1-4, 349-370. doi: 10.1016/0022-1694(95)92779-D.

TERZIC, J., STROJ, A. & FRANGEN, T. (2011): Hydrogeologic investigation of karst system properties by common use of diverse methods: a case study of Licka Jesenica springs in Dinaric karst of Croatia.- Hydrological processes. doi:10.1002/hyp.91

TOMASKO, D., FISHER, A., WILLIAMS, G. P. & PENTECOST, E.D. (2001): A statistical study of the hydrological character of the Edwards aquifer.- Argonne National Laboratory, University of Chicago for the U.S. Department of Energy, 34 p. doi:10.1061/40856(200)149.

VANDENSCHRICKA, G., VAN WESEMAELA, B., FROTA, E., PULIDO- BOSCH, A., MOLINAB, L., STIEVENARDC, M. & SOUCHEZD, R. (2002): Using stable isotope analysis (δD-δ18O) to characterise the regional hydrology of the Sierra de Gador, south east Spain. - J. Hydrol., 265, 43-55. doi: 10.1016/S0022-1694(02)00097-5.

VLAHOVIC, I., TISLJAR, J. & VELIC, I. (2009): Tercijarne karbonatne brece (paleogen-neogen- Pg, Ng) [Tertiary Carbonate Breccia, PGNg - in Croatian]. - In: VELIC, I. & VLAHOVIC, I. (2009): Tumac Geoloske karte Republike Hrvatske 1:300000 [Explanatory Notes of the Geological Map of the Republic of Croatia in 1:300.000 Scale]. Croatian geological survey, 141 p.

ZANINOVIC, K., SRNEC, L., PERCEC-TADIC, M. (2004): Digitalna godisnja temperaturna karta Hrvatske [A Digital Annual Temperature Map of Croatia - in Croatian].- Hrvatski meteoroloski casopis (Croatian Meteorological Journal), Zagreb, 39, 51-58.

ZUGAJ, R. (2000): Hidrologija [Hydrology - in Croatian].- Rudarskogeolosko- naftni fakultet, Sveuciliste u Zagrebu (Faculty of mining, geology and petroleum engineering, University of Zagreb), 407 p.

Manuscript received October 09, 2012

Revised manuscript accepted November 27, 2012

Available online February 28, 2013