Northern Velebit (Croatia) karst hydrological system: results of a preliminary 2H and 18O stable isotope study

Researching the hydrogeological properties of karst systems is very challenging due to their extreme heterogeneity. A grey-box approach in karst research combines the results from classical hydrogeological methods with direct observations within the karstic underground, i.e. in the caves. Isotope research has become a widely used method in the investigation of karst systems. The results presented are of a preliminary 2H and 18O stable isotope study of the Northern Velebit karst system (Croatia) employing the grey-box approach. Groundwater samples were collected during two summer expeditions in deep caves within the karst massif. Monthly precipitation samples were also collected (at three locations between approx. 900 m and 1600 m altitude), as well as water samples at some of the most significant springs, i.e. discharge outlets of the system. For a single expedition, the stable isotope composition is almost constant, i.e. the stable isotope measurements are within the measurement error across the complete cave profile. Simi­ lar characteristics across different caves during the same year were also noted. Samples of water from the springs were taken during base­flow conditions and they have similar isotopic contents to the cave water. The results obtained indicate that homogenization of the water already occurs within the subsurface epikarst zone above the sampling locations in the caves, but a future extended sampling campaign during variable hydrologic conditions is needed to confirm this. The final research goal is to establish a conceptual grey­box model for the functioning of this complex hydrogeological system. in karst hydrogeological research. It is based on rainfall – spring discharge modelling and water quality relationships, which are usually insufficient to reveal all the complexity of the various processes that control groundwater flow through the karst underground. The grey-box approach to karst systems research, besides the input-output relations, employs additional data acquired directly from underground i.e. in caves. Therefore, grey-box models are often more realistic and should be able to determine in more detail the system’s flow paths, linkages, stores, capacities, and throughput volumes than the black-box model (FORD & WILLIAMS, 2007). The research area of interest, the area of Northern Velebit (Croatia), has been the subject of earlier hydrogeological investigations (PAVIČIĆ, 1997; STROJ, 2010), but a number of parameters required to model the system’s behaviour are still unknown (MALOSZEWSKI, 2000). During the last few decades, many deep caves were discovered in the Northern Velebit karstic massif (BAKŠIĆ et al., 2013), providing an opportunity to gather some new and direct insights into its hydrogeology (STROJ & PAAR, 2019). The aim of the present preliminary study was to gather new insights into water infiltration and vadose flow properties within the karst massif by determining the δ2H and δ18O of water in various compartments of the system. Precipitation samples were collected from three locations situated on Mt. Velebit hilltops, while groundwater samples were collected in the caves and at the main discharge outlets of the system. The results presented provide guidelines for establishing a new research project supplemented Article history: Manuscript received March 01, 2019 Revised manuscript accepted June 10, 2019 Available online October 31, 2019


INTRODUCTION
Dissolution of soluble rocks such as limestone, dolostone and gypsum, forms a special type of hydrogeological system, which is often reflected at the surface by the distinctive landscape known as karst. Karst aquifers are characterized by extremely heterogeneous hydraulic properties manifested by rapid groundwater flows concentrated in conduits and disolutionally widened fractures, usually sparsely distributed, while the storage of by systematic and more frequent data collection. This could ultimately contribute to the development of a comprehensive conceptual hydrogeological model of the Northern Velebit system.

STUDY AREA
The study area is located in the Croatian part of the Dinaric karst, specifically in the northern part of a 145 km long Velebit coastal mountain range (Fig. 1), which is part of the external Dinaric mountains separating the Adriatic Sea coast from the large karst poljes of the Lika, and Gacka rivers that sink in its hinterland. The Northern Velebit area, designated a National Park, covers an area of approximately 109 km 2 . It has a form of high karst plateau Geological observations in caves revealed complex geological structures expressed by widespread occurrences of compact carbonate breccias down to depths of approx. 1000 m below surface, e.g. in the Lukina jama -Trojama and Velebita cave system (STROJ & VELIĆ, 2015). The breccias in the caves often irregularly intersect a normal succession of Jurassic limestones which points to their tectonic origin. This is commensurate with the assumption that the origin of the breccias observed at the surface is related to strong tectonic movements during the Upper Paleogene and Lower Neogene (BAHUN, 1974). Pronounced resistance to the mechanical erosion of massive carbonate breccia together with its solubility dominates the surface relief forming processes. This enhances the development of distinctive karst corrosion landscapes where the entrances of the most important caves are located. Such landscapes, together with occurrences of extremely deep caves, are instrumental in the National Park designation. Discovery and research of deep caves also revealed a very deep vadose zone of more than 1000 m inside the carbonate massif, as well as some characteristics of epiphreatic zone dynamics at the bottom of the vadose zone (BAKŠIĆ et al., 2013;STROJ & PAAR, 2019).
Sinkholes of the Lika and Gacka rivers are located on the eastern foothills of the Northern Velebit massif (Fig. 2). In the past, the rivers sank completely into the karst massif and reappeared at numerous coastal springs and vruljas (submarine (1) traced sinkholes in the massif`s hinterland, (2) main coastal discharge zones, (3) regional faults, (4) underground water flow connections determined by tracing tests, (5) sampled caves, (6) precipitation sampling locations, (7) meteorological station, (8) hydrogeological barrier formed by clastic rock formations, the rest of the area is built of various carbonate rocks (adapted with permission from STROJ, 2010).
intersected by numerous karst depressions separated by peaks and ridges (the highest peaks are 1600 -1700 m a.s.l.). The largest part of Northern Velebit is composed of Jurassic carbonate rocks, mainly limestones, and massive calcareous (limestone) breccias of Upper Paleogene to Lower Neogene age (MAMUŽIĆ et al., 1969;VELIĆ et al., 1974;VELIĆ & VELIĆ, 2009). Clastic rock formations of Triassic and Palaeozoic ages form a hydrogeological barrier in the middle-southern parts of the Velebit mountain range. The deeper position of these rocks in the northern part enables deep karstification and underground water flow from the hinterland towards the Adriatic Sea (Fig. 2).
The geological structure of the Velebit mountain range is the consequence of two principal periods of tectonic activity (BLAŠKOVIĆ, 1998;KORBAR, 2009). The first was the main phase of the Dinarides orogenesis (Eocene -Miocene) with compressive movements oriented NE -SW. As a result of the regional tangential stress, deep nappe structures, folds and regional faults of Dinaric strike (NW -SE) have been formed. The second, socalled neotectonic phase, started probably during the late Miocene or early Pliocene, when the main stress changed to NS, resulting in further uplifts and transpressive deformations of older structures. The later phase played a crucial role for the present mountain morphology. springs) along the coast on the western side of the massif (Fig. 2) (PAVIČIĆ, 1997;STROJ, 2010). The distance from the sinkholes to the Adriatic Sea is approx. 20 km. Today, the water from the two rivers is used for electricity production and is routed via a tunnel to the power plant on the coast. Only occasional water surpluses sink into sinkholes occurring on average only a few days per year. Taking into account the present conditions, the precipitation in the massif area is assumed to be the main water input for the system.

SAMPLING AND DATA
Precipitation samples were collected on a monthly basis from May 2012 to September 2013 on Velebit at three locations -Oltari, Sića, and Crikvena (Fig. 2). Geographical coordinates and altitudes of the precipitation sampling sites are provided in Table 1. Precipitation sampling was performed usingsamples with antievaporation system (GRÖNING et al., 2012). A total of 35 monthly precipitation samples were analysed for their stable isotope content ( Table 2). Some of the samples are missing due to  sampling locations inaccessible during severe winter weather conditions. Also, some samples had to be discarded due to an overflow of collectors (November 2012), while others were not suitable for the analysis due to insufficient volume (August 2012, summer 2013). Water samples in the Velebita cave system (VCS, entrance approx. 1550 m a.s.l., Fig. 2) were collected during Summer 2012 and Summer 2013 spelaeological expeditions. A total of 20 groundwater samples were collected in the caves. Samples were collected during dry summer days as summer storms increase the probability of sudden, increased water flows through cave passages, which is dangerous for cavers. The VCS sampling in 2012 started at a depth of approx. 100 metres, and continued to 1026 metres depth (Table 3). In 2013, two samples were collected in the VCS at a depth of approx. 600 m, and one sample was collected in the Lukina jama cave (entrance at approx. 1450 m a.s.l., The main coastal springs along the western foothills of the massif: Bačvica (BC), Jurjevska Žrnovnica (JZ) and Starigrad Senj (SS), were sampled in the summer of 2012 (Fig. 2). Additionally, two samples were collected from Štirovača spring (ST), located in the central part of the Velebit massif at 1100 m a.s.l. (Fig.  2). A total of 9 spring water samples were collected (Table 4).
All collected samples (precipitation, groundwater collected in caves and spring water) were stored in 50 ml, double capped, polyethylene bottles until analysis.
The Meteorological and Hydrological Service provided meteorological data for the Zavižan meteorological station (1600 m a.s.l.), the only permanent meteorological station on Northern Velebit.

MEASUREMENT AND ANALYSIS
Stable isotope measurements were performed using the water equilibration technique (HORITA & KENDALL, 2004) on a Delta plusXP (Thermo Finnigan) isotope ratio mass spectrometer. An HDO eq48/24 (IsoCal) equilibration unit and a Dual Inlet   CRISS et al., 2007). Normalization and analyses of measurements were performed by the US Geological Survey Laboratory Information Management System. Statistical calculations were conducted in Statistica 13.0 (Statsoft Inc., Tulsa, OK, USA). The descriptive statistics for data that are normally distributed are presented with arithmetic mean and standard deviation, while data that do not follow a normal distribution are presented with the median and minimum and maximum values. To test whether a sample distribution matches the characteristics of a normal distribution we used the Kolmogorov-Smirnov test. A measure of the linear correlation between two variables is expressed with the Pearson correlation coefficient (R). The relationship between two variables was determined by linear regression analysis. Testing the differences between the continuous variables was accomplished by t-test (for normally distributed data) or Mann-Whitney test (for data that do not follow a normal distribution), and regarding the categorical variables by c 2 -test. The level of statistical significance was chosen at p < 0.05. do not differ significantly (Mann-Whitney test, p = 0.43). The average monthly air temperatures for the period showed typical seasonal behaviour with the highest temperatures in summer and the lowest temperatures in the winter months (Fig. 3A).

RESULTS
Monthly precipitation during the hydrological year 2013 (median 218.6 mm; minimum 45.1 mm; maximum 462.9 mm) was significantly higher compared to the hydrological year 2012 (median 86.2; minimum 0.8 mm; maximum 278.3 mm) (Mann-Whitney test, p = 0.04). As shown in Fig. 3A, this is particularly evi-dent during the period Jan -Mar 2012 (very low precipitation) and Jan -Mar 2013 (rainy period).
During the hydrological year 2013 snow was present for longer in comparison to the previous hydrological year (c 2 -test, p = 0.003, Fig. 3B). In the hydrological year 2013 the snow cover was continuously present for 164 days (melting in May 2013), while during the hydrological year 2012 the snow cover was present for 100 consecutive days (

Precipitation
Due to the small number of precipitation samples, we weren't able to set up a Local Meteoric Water Line (LMWL) that would be relevant for a longer period. Nevertheless, there is a statistically significant correlation between δ 2 H and δ 18 O precipitation values (R = 0.99; p < 0.001). A corresponding precipitation regression line: d 2 H = 7.15 • d 18 O + 6.14 ‰ (R 2 = 0.98) plots close to the Global Meteoric Water Line (GMWL) that expresses the worldwide average (CRAIG, 1961) (Fig. 4).
Given the strong correlation between δ 2 H and δ 18 O precipitation values, only δ 18 O values are mentioned below. The time series of precipitation δ 18 O shows higher values in the summer months than during the rest of the year (Fig. 5A). This is typical for low latitudes where 'the amount effect' accounts for δ variations: low isotopic signals in months with significant precipitation, and more positive isotopic signals in months with scarce rainfall (DANSGAARD, 1964).
The average d-excess value of precipitation samples is (13.16 + 2.72) ‰. Fig. 5B shows that precipitation d-excess values are lower during summer than during the rest of the year. This is comparable to previously reported results from the bordering region with similar climatic conditions . Therefore, we can conclude that the sources of the precipitation are comparable, i.e. that the dominant summer precipitation in the area of Northern Velebit originates in the Atlantic Ocean, while in winter there is the predominance of Mediterranean-originated precipitation.

Groundwater
The isotope measurements for samples collected in 2012 across the VCS profile (depths >1000 m) show a very small variation of δ 18 O values ( Figure 6, Table 3). The average δ 18 O value of these samples is (-9.59 + 0.13) ‰. The values are within a narrow range (< 0.5 ‰) and without the presence of a uniform vertical gradient. Stable isotope composition of the samples collected in caves corresponds to the lower part of the precipitation cluster (Fig. 4), indicating that the groundwater flowing through the cave during low water conditions is dominantly fed by the September -May precipitation (Fig. 5).
The relatively constant isotopic composition of the samples suggests that groundwater homogenization already occurs within the epikarst zone. The epikarst zone comprises a relatively shallow, but most intensely fractured and karstified subsurface zone where precipitation quickly infiltrates and could be temporarily stored (MANGIN, 1974;WILLIAMS, 1983). Water stored in such epikarst storages during the wet season slowly drains into the underlying vadose zone during the dry summer season. Therefore, we can conclude that consistent water characteristics across the VCS profile suggest that the cave streams are dominantly fed from the epikarst storages. The substantial inflow of water with different characteristics from deeper storages, i.e. from fractures in rock walls along the cave profile, would probably be reflected by more pronounced variations in the stable isotope composition with cave depth. It should be mentioned once again that all samples were taken during low water conditions. If samples were collected during summer storms, when fast infiltration can bypass epikarst storages through fractures widened by dissolution, δ 18 O would probably show different values, i.e. the isotopic composition would be more influenced by recent (summer) precipitation.
Cave samples collected during the 2013 expedition have lower δ 18 O values than those from the 2012 expedition (Fig. 6,  Table 3). These observed differences are tentatively explained by the characteristics of the rainfall period that preceded the second expedition. Namely, lower isotopic values could be associated with higher levels of snow and rainfall during the hydrological year 2013 in comparison with hydrological year 2012 (Fig. 3).
Unfortunately, sampling along the complete VCS profile was not performed in 2013, and only samples at approx. 600 m depth were available. It is interesting that the sample taken in the Lukina jama in the same year has a very similar isotopic value. This suggests similar recharge mechanisms for different caves in the area.
It is generally considered that the preserved seasonal amplitude of the spring water isotopic signal indicates a mean residence time of < 1 year (TRČEK & ZOJER, 2009, LAUBER & GOLD-SCHEIDER, 2014. If lower isotopic values of the sampled groundwater in 2013 (Fig. 6) are a consequence of the characteristics of the preceding wet season, this points to a relatively short mean residence time of groundwater in storages before draining to the cave conduits. Nevertheless, this should be confirmed by longer and more extensive sampling of the precipitation, and groundwater in the caves.
The sample collected at a depth of 45 m from the Sirena pit in 2013 was depleted in comparison to other samples collected   that year (Fig. 6). We may attribute the depletion of the sample to the probable origin of the water sampled from the snowmelt, as for snow lower isotopic values are expected (MANCE, 2014). Also, perennial snow accumulation was reported by the cavers in a deep dolina where the cave entrance is located.
For d-excess values of the samples collected in caves, we may also see that most of the data plots within the range of the measurement error (Fig. 7). The average d-excess value of cave collected samples is ( (t-test, p = 0.28). This is an additional indicator of predominant groundwater recharge by the winter precipitation. The extreme d-excess values have samples collected at the highest altitudes i.e. the lowest cave depths. This may be explained by some contribution of newly infiltrated precipitation in the sampled water. Fig. 4 shows that the samples from springs plot above the samples collected in caves. Water from coastal springs was slightly brackish due to seawater intrusion into the karst aquifer, which was manifested by the elevated conductivity of the spring water. Specific conductivity measured during sampling at the BC, JZ and SS springs was 6.7-7.0 mS/cm, 7.9 mS/cm and 13.8-14.7 mS/cm, respectively. Specific conductivity of the spring water during high flow conditions (when seawater intrusion doesn't affect it) was approx. 0.3-0.4 mS/cm (STROJ, 2010). We estimated the seawater proportion in the spring water samples by applying two-component mixing analysis (PETERLLA & CELICO, 2013) based on the sample conductivities, unaffected groundwater conductivity and seawater conductivity (55 mS/cm, measured during sampling). The calculated proportion of seawater in the spring water was approx. 12 % at BC, 14 % at JZ and 25 % at SS springs. We further applied two-component mixing analysis for corrections of coastal springs' isotopic compositions. Corrections were based on the estimated proportion of seawater and presumed seawater isotopic composition of δ 18 O ≈ +1‰ (GAT et al., 2003). After the correction, groundwater samples from the springs have similar values to water from the caves in the same year (Tables 3  & 4), confirming the coastal springs are predominantly fed by the recharge from the massif.
Stable isotope compositions of groundwater collected in both, caves and springs, coincide with the precipitation regression line (Fig. 4). This indicates the meteoric origin of the groundwater without significant influence of additional post-precipitation fractionation.
The only sampled inland spring, located in the central, but somewhat lower part of the massif (1100-1300 m a.s.l.) is the Štirovača (ST) spring (Fig. 2). Both ST samples have δ 18 O values which are comparable to, but slightly less depleted than those collected in caves (Tables 3 & 4) suggesting it is fed from its immediate surroundings. This is consistent with the geological situation: the spring is located close to the geological boundary with the impermeable clastic rocks situated on the southern margin of the Northern Velebit.
The δ 18 O values of samples collected in the caves and the ST spring are comparable to values for some of the Gacka river springs sampled during summer months of 2005(MANDIĆ et al., 2008OZYURT et al., 2014). This indicates similar mechanisms of winter (wet) season recharge of water storages, which are slowly drained during the dry season. The ST d-excess values are comparable to d-excess values of samples from the caves (Tables 3 & 4). This is another indication of similar water feeding into the caves and the ST spring. Lower d-excess values of coastal springs are probably a consequence of mixing with seawater or/and a small contribution from recent summer precipitation.

Preliminary conceptual model of the Northern Velebit hydrogeological system
Taking groundwater samples directly from the deep caves within the massif permitted the use of the grey-box model instead of the black-box model that takes into account only system's input and output data. Fig. 8 shows a conceptual model of the Northern Velebit hydrological system developed on the basis of the results presented here, including the results of previous spelaeological and hydrogeological research (PAAR et al., 2008;STROJ, 2010;BAKŠIĆ et al., 2013;STROJ & PAAR, 2019). The proposed conceptual model includes provision of significant storage capacities within the epikarst zone, situated immediately below the surface, down to the depth of a few metres (or few tens of metres in case of massive and well lithified carbonate breccias). During the dry season, the flow in cave conduits is mostly fed by the water stored in the epikarst during the previous winter / wet season. Across the vadose zone, the flow is fast and predominantly concentrated in conduits, which is expressed by constant isotope characteristics. The preliminary findings presented earlier suggest that the mean residence time of water stored in the system is relatively short, possibly less than a year. This is indicated by the cave water stable isotope composition influenced by characteristics of the preceding winter / wet season. Coastal karst springs, that are the main discharge points of the system, are predominantly fed by the water of similar characteristics as that observed in the caves. The infrequent nature and short duration of concentrated inflows from sinkholes on the Lika and Gacka rivers under present conditions limits their significance. Therefore, it may be concluded that epikarst storages on the massif are probably the main source of spring water during the base flow conditions.

CONCLUSIONS
The paper presents the results of a preliminary stable isotope study of the waters of the Northern Velebit karst hydrological system. The aim of the study was the application of a grey-box modelling approach to karst hydrogeological research. The approach combined the results from classical hydrogeological methods with the findings from the sampling of the water in the caves.
The study used stable isotope analysis of water samples as a tool to investigate the functioning of the karst hydrogeological system.
The results of this study, together with the results from previous research enabled us to develop a preliminary conceptual recharge and storage model of the Northern Velebit karst system. The main features incorporated in the model are: • groundwater storage is mainly related to the function of the epikarst zone; • recharge of the epikarst occurs mainly during the wet / winter season; • the mean residence time of epikarst storage is relatively short, probably within one year; • the deeper vadose zone has a low storage capacity and water flows within are concentrated in conduits; • system outflow via coastal springs and vruljas (submarine springs) in the base flow conditions is dominantly fed by epikarst storage of the massif. Based on the obtained results, we conclude that the grey-box modelling approach is an improvement in the investigation of complex karst systems. In addition to classical hydrological methods, spelaeological investigations play an important role in this type of approach and give researchers the opportunity to study different hydrological processes in situ.
Although the number of analysed water samples is relatively low, limiting their interpretation, this study provides preliminary results, enabling the establishment of a recharge and storage model for the system that will serve as a basis for future research.
Establishing a future longer term monitoring programme for the Northern Velebit hydrogeological system would provide data required to further test and upgrade the presented model. This would provide a better and more detailed understanding of the functioning of the Northern Velebit karst system and enable efficient water management and nature conservation to be undertaken.