Hydrological Behavior of Karst Systems Identified by Statistical Analyses of Stable Isotope Monitoring Results

The article presents findings of a two-year systematic study of stable isotope content in two karst groundwater resources in Primorsko-goranska county (Croatia): the Martinšćica wells (MWs) and the Dobrica spring (DBC). The temporal and spatial variation of hydrogen and oxygen isotopes is commonly studied in conjunction with hydrogeological conditions such as groundwater dynamics and discharge conditions. However, since this information was incomplete, we were forced to work with limited data and rely on analyses of stable isotope monitoring results. The obtained results show that winter precipitation is the most common recharge source for the systems, and the average residence time of water in the subsurface is less than a year. Furthermore, the MWs system is a typical dual-porosity system with dominant base flow. The results of the nonparametric regression analysis show that the possibility of seawater intrusion into the spring affecting DBC isotope content cannot be ruled out. We believe that the results presented in the paper demonstrate that when combined with statistical analyses, environmental stable isotopes are a powerful tool for gaining insights in


Introduction
It is estimated that 20-25% of the world population's drinking water comes from karst systems [1][2][3]. Groundwater resources from karst aquifers, in particular, are a vital natural resource in the Mediterranean region, as karst springs constitute many Mediterranean countries' only supply of water [4,5]. Croatia is no exception, as settlements in Croatian karst, which cover half of the country's land, are almost entirely supplied with drinking water from karst aquifers.
Knowledge of spatial and temporal variations in the functioning of karst aquifers is the basis for sustainable water management, allocation of water resources, and their protection. In Croatia, as in most Mediterranean countries, this is especially true during the dry summer months, when the population on the coast increases greatly due to tourism. During this time, anthropogenic pressure on karst aquifers increases due to both possible overexploitation and an increased risk of pollution [6][7][8][9].
The study of karst systems is challenging because they are characterized by heterogeneous physical, chemical, and biological properties, and complex groundwater flows [10]. In addition, karst areas are subject to constant change due to dissolving water activities, expansion of fissure systems, soil erosion, and collapses of underground cavities [11]. Due to their complexity, interdisciplinary research on the hydrology of karst aquifers has been practiced and proposed as a model for future directions in hydrology [12,13].
The use of stable environmental isotopes to complement conventional hydrogeological methods contributes to a better understanding of karst hydrogeology [14,15]. The stable water isotopes 18 O and 2 H, which are components of the water molecule, diffusely percolate/seep throughout the catchment during rainfall and thus can be used as ideal environmental traces for studies of karst aquifers [16][17][18].
Our goal was to characterize the functioning and origin of groundwater from two adjacent unconfined karstic aquifers in Kvarner Bay (Croatia). We had to rely on stable isotope analysis for our research because there was limited available hydrochemical data and no appropriate data on flows/dynamics of water discharge and groundwater levels for the analyzed water sources.
In addition to gaining a better understanding of the behavior of karst aquifers under various hydrologic conditions, the objective of the study was to evaluate the extent to which statistical analysis of stable isotope content in water is consistent with conclusions drawn from a traditional groundwater physicochemical parameters. Confirmation of good agreement would indicate that thorough statistical analysis of δ 18 O and δ 2 H time series is suitable for determining the hydrological behavior of karst systems when conventional parameters are not available.

Study Area
The study area is located in the Dinaric Karst region in western Croatia ( Figure 1, upper right corner). It includes three zones of groundwater discharge in the wider area of the Gulf of Rijeka: the spring zone within the city of Rijeka along the right bank of the Rječina riverbed, the Martinšćica Bay area, and the spring outflow zone in the Bakar Bay. These complex and interconnected recharge zones are generally treated as the unified Rijeka-Bakar groundwater zone [19]. The area of this zone is about 600 km 2 . It is highly fractured and extends over the area between the sea and the mountain massif with peaks over 1500 m above sea level. The highest parts of the basin are also characterized by the highest average annual precipitation in Croatia, up to 3500 mm. The average annual precipitation for the entire groundwater body is about 2170 mm, the average annual air temperature is 8.7 • C, and the average annual runoff is 30.5 m 3 s −1 [19]. The geological structure of the area is complex. The numerous lithological components of the area can be divided into three basic groups: carbonate rocks, flysch, and Quaternary clastic deposits. Most of the area is composed of carbonate rocks, whose fissure systems range from millimeter systems to channel and cave systems. The flysch deposits are of low permeability and are a barrier to groundwater flow. Clastic deposits have low water retention capacity and their impact on groundwater dynamics is not significant (Figure 1).
The main distribution of groundwater from the mountainous hinterland to the springs in the coastal area takes place on the plateau of Grobničko polje at an altitude of 300 m above sea level. On the northern edge there are a number of occasionally active springs ( Figure 1). On the southern edge, there are swallow holes (ponors) from which groundwater flows to the coastal springs of Zvir (ZV), the springs in the Martinšćica Bay, and the springs in the Bakar Bay [21,22].
For most of the year, the drinking water supply of the city of Rijeka and surrounding settlements is ensured by the largest spring in the region: the Rječina Spring (RJ). RJ is an overflow spring that regularly dries up during the summer. During this period, the local public utility responsible for the water supply, has to provide sufficient water not only for the city of Rijeka, but also for the surrounding settlements, which are all popular tourist destinations. The water is then diverted from the largest permanent spring in the area: ZV. The same public utility manages Martinšćica wells (MW). Since the early twentieth century, approximately ten wells have been excavated in the Martinšćica valley, which is located at an elevation of 3-5 m above sea level. The depth of the wells is approximately 6 m, and only five of them are currently operational ( Figure 1). MWs are in an area where, after heavy rains, there are occasional significant leaks of groundwater from the marginal parts of the same basin that feeds both RJ and ZV.
The public utility also manages the water supply from three nearby aquifer springs: Perilo (PER), Dobra (DB), and Dobrica (DBC). These springs are located in Bakar Bay very close to sea level. MWs and springs in Bakar Bay piqued our interest due to their location. MWs are found in a marginal and underexplored area of the Rječina basin, next to the Bakar springs basin, which is linked to the sinking zones of Grobničko polje. In this article, we will look at one of the MWs (specifically, well 2-W2) from the Rječina aquifer and DBC as a representative of the Bakar Bay springs.
The amounts of surface and groundwater flowing into the sea are registered in the Martinšćica hydrological station (MAR, Figure 1). It is a hydrological system with significantly lower flow rates than at the main groundwater collector of the Rječina River. The average annual flow at the RJ is 7.4 m 3 s −1 , while the flow for MAR is about 1.3 m 3 s −1 [23]. In the event of heavy precipitation, the flow recorded on MAR includes both groundwater discharging at the location and surface waters from the Javor stream ( Figure 1). The waters of Javor are collected in its flysch basin beneath which flows groundwater that recharges the Martinšćica spring zone. Because these are marginal parts of aquifers near the sea coast, the range of seasonal groundwater level fluctuations in the outflow zone in the Martinšćica area is up to about 5 m, whereas the range of seasonal groundwater level fluctuations in the deeper hinterland, in the Grobničko polje area, is up to about 50 m. Unfortunately, no discharge measurements were taken at the DBC. The only information we have is the average annual flow of the source runoff zone in Bakar Bay, which is approximately 0.34 m 3 s −1 [19]. The dynamics of groundwater fluctuations and discharges in the Grobničko polje area are linked to groundwater in the Bakar area, which has a small range of groundwater level oscillations (order of 1 m) and discharges directly into the sea [23].
Problems with salinization of spring waters occur at the southeastern coastal margin of the Rječina River aquifer. In the Bakar Bay area, salinization is caused not only by the overexploitation of water reserves, but also by the dynamics of groundwater discharge in contact with the sea [24].

Sampling and Data
Groundwater samples were collected weekly for two consecutive years (April 2010-April 2012). During the first sampling period, 55 samples were collected at each of the five MWs. The analysis showed that there was no significant difference in the abundance of isotopes 18 O and 2 H in the water of the sampled wells [25]. Therefore, only W2 was included in further sampling because it was considered by the water utility company to be the most important for water supply. During second sampling period, an additional 33 samples were collected from W2. A total of 100 samples were collected from the DBC spring. Samples were stored in HDPE bottles with double caps until analysis. To gain a better understanding of the hydrologic behavior of W2 and DBC in the analysis presented here, previously published data for RJ were also used [17].
To compare the isotopic content of precipitation and groundwater, we used a set of precipitation isotope data from the Kukuljanovo (KUK), Škalnica (SKAL), and Platak (PLAT) rain gauging stations ( Figure 1). Part of these data has already been used for hydrogeological analysis of RJ, ZV, PER, and DB [17]. Geographical coordinates and elevations of precipitation sampling sites are given in Table 1. Most of the presumed recharge area includes an inaccessible mountainous area where no additional precipitation gauges could be installed. Cumulative monthly samples were collected in 3.5-L rain gauges to which 100 mL of paraffin oil was added to prevent evaporation. After separating the rainwater from the oil, the rainwater samples were stored in the same manner as the groundwater samples. Isotopic composition was measured by the water equilibration method using an isotope ratio mass spectrometer (IRMS, DeltaplusXP, Thermo Finnigan, Bremen, Germany) in conjunction with an equilibration unit (HDOeq, IsoCal). The results are expressed according to the VSMOW2-SLAP2 scale. The measurement accuracy achieved was better than 0.1‰ for δ 18 O and better than 1‰ for δ 2 H. Groundwater temperature was measured in situ, while measurements of water turbidity, pH, electrical conductivity (EC), and Cl concentration were carried out in the laboratory of the municipal water utility in Rijeka.

Meteorological and Hydrological Situation
To better understand the changes in the isotopic composition of groundwater, we examined the meteorological and hydrological situation during the study period ( Figure 2). For this analysis, we used data from the official meteorological station of the Croatian Meteorological and Hydrological Service (CMHS) in Kukuljanovo. We chose this meteorological station because it is located in the middle of the basin and is closest to our rain gauge station with the longest sampling period (KUK). Martinšćica hydrological station (MAR, Figure 1) provided data for the analysis of the hydrological situation. It should be noted that the discharge data for years 2007, 2008, 2017, and 2018 are incomplete.  It should be noted that the sampling period included two extreme years in terms of precipitation. There were three months in 2010 with precipitation above the multi-year average: September, November, and December. This rainy period was followed by a period with precipitation that was significantly lower than the multi-year averages (Kruskal-Wallis test, p = 0.019). Precipitation was completely absent in August 2011 and March 2012.
MAR flows were higher than ten-year averages in the fall and early winter of the first year of sampling (September and December 2010 and January 2011, Figure 2). The average MAR discharge in the second year of sampling was not only statistically significantly lower than the previous year, but also lower than ten-year averages (Dunn test, p < 0.016).

Statistical Analysis
We used the auto-correlation function (ACF) to analyze the mean daily MAR discharge. ACF is a collection of auto-correlation coefficients r(k) xx that vary with time lag k: n is the number of data points for which m autocorrelation coefficients are computed. Mangin (1984) states that m should not be greater than n/3 [26]. x t = (x 1 , . . . , x n ) is a time series for which auto-correlation coefficients are calculated, and x and σ xx denote the mean and standard deviation for that time series.
The shape of the ACF discharge graph can be used to infer the karst aquifer's retention capacity and degree of karstification. The time required for r(k) xx to fall below 0.2 is known as the karst system's memory effect [26]. The karst aquifer with a longer memory effect has greater retention capability and a less developed karst aquifer conduit system than the karst aquifer with a shorter memory effect.
For statistical analysis of groundwater stable isotope content, we used Gaussian Mixture Modeling (GMM) in addition to standard empirical methods such as descriptive statistics, correlation analysis, and the t-test. GMM was used to find the underlying densities of time series based on means, standard deviations, and weighting factors. GMM has been successfully used in the analysis of stable isotopic composition of karst spring water [25]. The basis of GMM is the assumption that the cumulative probability density function g(x) of the isotope series is a weighted sum of k underlying Gaussian components: where i = 1, . . . , k are components, i.e., different water masses in the system. f i is the probability density of the normal distribution corresponding to a single water mass. Accordingly, µ i is the mean isotopic value of the respective water mass and σ i is the corresponding standard deviation. The weighting factor w i indicates the proportion of the contribution of each component to the total isotopic composition. Bayesian information criterion was used to determine the number of underlying components (k). The analysis was performed using the R package mclust [27].
To determine the relationship between DBC water δ 18 O levels and multiple predictors, we used multivariate adaptive regression splines (MARS). MARS is a nonparametric regression technique that approximates the relationship between the dependent variable and multiple parameters through piecewise regression [28]. MARS predicts a function by linear combinations and interactions of adaptive piecewise linear regression, i.e., the "basis function (BF)". Therefore, f (x) (in our case δ 18 O) of the model MARS can be specified as follows [29]: where x is the independent variable, β 0 is a constant, λ i (x) is BF, β i is the coefficient of the i th BF, and n is the number of BFs in the model. The least squares method was used to calculate the coefficients. The BFs are functions of the form max(0, x − α), where α is constant corresponding to a knot. To obtain the BF, two adjacent splines intersect at a knot.
The stepwise forward method and the stepwise backward algorithm were used to build the MARS model. The stepwise forward approach adds BFs to Equation (4) and searches for potential knots to improve model performance. However, getting too many BFs with this method can cause the model MARS to overfit. The stepwise backward method is used in the second phase to address this problem. In this phase, redundant BFs with the smallest contributions are removed from the BFs used in the stepwise forward approach to determine the best submodel. Generalized cross-validation (GCV) is used to eliminate redundant BFs from the MARS model. It is calculated as follows: where N is the sample size, y i is the observed value of the response variable ( is the MARS predicted value, M is the number of BF, and d is the penalty parameter. It is considered that the optimal d value is in the range of 2 to 4 [29]. MARS is increasingly used in environmental studies because the relationships between environmental parameters are often non-linear [30,31]. MARS analysis was performed in the R package earth [32].

Discharge Analysis
As previously stated, no water level or discharge measurements are taken at the MW2 or DBC. The only data available are discharges recorded at the MAR hydrological station, which include both groundwater and surface water flows. To gain insight into the water dynamics, we compared mean monthly discharges on RJ (as a representative of the upper aquifer part) and MAR from 1975 to 2016. Due to a lack of data, the years 2001, 2007, and 2008 were excluded from the analysis.
MAR's mean monthly discharge is lower than RJ's (Mann-Whitney, p = 0.0007). The discharge graph shows a seasonal pattern, with summer minimums and fall maximums. MAR lacks a distinct discharge peak in April, as seen in RJ, due to its lower elevation and thus lower snowmelt component in the discharge (Figure 3). The analysis of autocorrelation functions of mean daily discharge (Figure 4) shows that the MAR system has a longer memory effect (39 days) than RJ (22 days). As a result, the MAR system has a more prominent base-flow component. Drying of the groundwater sources caused the occurrence of negative autocorrelation values, with RJ drying out more frequently than MAR.   Table 2). A common procedure in isotope hydrology is to calculate the deuterium excess (d-excess = δ 2 H−8·δ 18 O, [33]). The d-excess value of a single atmospheric precipitation event is mainly related to the temperature of the vapor-causing sea surface, i.e., to the evaporative conditions of the source of the precipitation. As we can see in Figure 6, the dexcess in precipitation shows a distinct seasonal behavior with lower values in summer and higher values in winter. Low precipitation events and events occurring during relatively hot weather conditions have low d-excess values [34], so this seasonality was to be expected. To avoid the influence of such low-precipitation events, we calculated amount-weighted d-excess values. The weighted mean d-excess value in the period coinciding with the time of groundwater sampling is 12.9‰ for KUK, 14.71‰ for SKAL, and 14.66‰ for PLAT. Means, ranges, and standard deviations of groundwater isotopic contents are shown in Table 3. Data for RJ were added for comparison [25]. The δ 18 O values of W2 water varied from −9.52‰ to −6.63‰ and the δ 2 H values ranged from −61.4‰ to −39.7‰ (Table 3). During the same period, the δ 18 O values of the DBC spring water varied from −9.2‰ to −7.1‰, while the δ 2 H values ranged from −53.5‰ to −42.0‰ (Table 3).

Stable Isotope Content of Precipitation and Groundwater
The linear relationship between δ 2 H and δ 18 O in natural waters was first recognized by Craig [35] and the corresponding regression line was named Global Meteoric Water Line (δ 2 H = 8·δ 18 O + 10‰, GMWL). In isotope hydrology, it is common to define this relationship at the local scale and construct corresponding local meteoric water lines, i.e., LMWLs. Previous analyses from the research area show that the isotopic content of the local groundwater plots is above corresponding LMWL [17].  If we compare the isotopic values of groundwater (Table 3) and monthly precipitation ( Table 2), we can see that the average values of groundwater are more negative than the average values of precipitation. We constructed two LMWLs to test the idea that winter precipitation primarily recharges groundwater: one for the warm months (April-September) and one for the cold months (October-March) of the hydrological year. For this purpose, we used the isotopic values of KUK precipitation (Table 4), since the longest measurements are available for this station.  Figure 7, we see that the isotopic values of the groundwater agree well with the LMWL for the cold season, confirming our assumption that the groundwater is fed predominantly by winter precipitation. Local groundwater lines (LGWL) for MWs (δ 2 H = 8.12·δ 18

Analysis of δ 18 O-Groundwater Time Series and Isolated Hydrological Events
Because of the high and significant correlation (R 2 > 0.93; p < 0.001) between δ 2 H and δ 18 O-groundwater values, for further analysis we focused only on δ 18 O. To gain a better insight into the situation, we included the RJ in the analysis as the highest yielding source in the area. Although the discharge of MAR is significantly lower than that of RJ, we can see that the leaps in the discharges occur at approximately same time (Figure 8).
Regarding the isotopic composition, the isotopic values of W2 and RJ are comparable in the summer months and in times of drought (Figure 8), when the base flow at the springs is drained and the water reserves are depleted. Therefore, we cannot reject the hypothesis that RJ and W2 share common water storage from deeper parts of their aquifer in the mountainous hinterland.
In Figure 8, we see that the shift toward less negative δ 18 O values coincide with heavy precipitation, generally occurring earlier in RJ than in W2 and DBC, and also that these less negative values remained longer in W2 and DBC than in RJ. The faster RJ reaction is most likely due to this spring's higher hypsometric position (325 m a.s.l.) and more pronounced fast-flow discharge component. This is consistent with the ACF analysis, which revealed that the lower part of the basin has a longer memory effect than RJ, i.e., a more pronounced base-flow component. Water on the RJ accumulates, rises, and flows at the flysch barrier's very edge, whereas groundwater flowing toward the MWs and DBC flows below the barrier that slows it down. In addition, some of the groundwater feeding the MWs and DBC comes from precipitation infiltrated into the karst aquifer in the lower parts of the basin. From the analysis, we conclude that shifts toward less negative δ-values are associated with the presence of newly infiltrated precipitation in the MWs system, while in the case of DBC they could also be caused by seawater intrusion into the spring [24].  (Figures 8-10).
In the event E1, a large amount of rain fell in a short time after a long period of drought. The δ 18 O-rainfall values for September 2010 are: −5.65‰ for KUK, −5.25‰ SKAL, and −5.35‰ PLAT. This resulted in changes in groundwater flux and isotope values, but also in large changes in the characteristics of several other water quality indicators (Figures 9 and 10). The increase in the proportion of newly infiltrated rainwater in the groundwater caused the groundwater temperature to rise. A sharp increase in turbidity was also observed in the studied locations. The reason for this, in addition to the presence of a proportion of newly infiltrated water, is the intensification of groundwater flow, i.e., its velocity and thus kinetic energy. This caused the movement and flow of previously deposited suspended solids in the karst aquifer [36].
In event E1, the newly infiltrated water caused a decrease in groundwater pH, but also an increase in EC. This suggests that the rise in EC in September 2010 may have been influenced by infiltration of more chemically aggressive soil water following intense rainfall, but also by flooding of previously drained channels [37]. This model of aquifer behavior, i.e., the described changes in groundwater properties, was not repeated in the same way in the other isolated situations (E2-E4), probably because the change in hydrologic-hydraulic conditions was not as intense.  To confirm the relationships observed in events E1-E4, we performed a correlation analysis. Tables 5 and 6 show the correlation coefficients for the groundwater properties studied at W2 and DBC, respectively. Table 5. Matrix of linear correlation coefficients between groundwater parameters of Martinšćica well 2 (W2). Only statistically significant coefficients are given. t w -water temperature, EC-electrical conductivity, NTU-turbidity; Cl-chloride ion concentration; t a -average air temperature between two consecutive samplings; mm-total rainfall between two consecutive samplings.  Table 6. Matrix of linear correlation coefficients between groundwater parameters of Dobrica spring (DBC). Only statistically significant coefficients are given. t w -water temperature, EC-electrical conductivity, NTU-turbidity; Cl-chloride ion concentration; t a -average air temperature between two consecutive samplings; mm-total rainfall between two consecutive samplings. As for the W2 (Table 5), there is an almost perfect correlation between the δ 2 H and δ 18 O values. Moreover, a statistically significant negative correlation was found between isotopic composition and air temperature, confirming that summer groundwater (which is recognized as the base flow) has the most negative values. Although a higher EC usually indicates "old" water with more solutes, in our case it was found that there is a statistically significant positive correlation between conductivity and turbidity, which we have already explained in the previous section by infiltration of chemically aggressive soil water following intense rain and flooding of previously drained karst channels. Turbidity is positively correlated with both precipitation and δ 18 O. Since in most cases the δ 18 O of the rainfall had less negative values than that of the groundwater, the suggestion that the shift in the isotopic composition of the groundwater toward less negative values is a result of new water infiltrating the system cannot be dismissed.
The correlation in Table 6 shows that we also have an excellent relationship between δ 2 H and δ 18 O in the case of DBC. There is a significant positive correlation between EC and Cl in DBC, which corresponds to marine intrusion into the spring. Isotopic composition is negatively correlated with air temperature, as in the case of W2, again suggesting that the base water at the source is associated with more negative values of δ 2 H and δ 18 O.

Gaussian-Mixture Modelling (GMM) of Groundwater Stable Isotope Data
GMM, applied to the frequency distributions of EC, has been used for some time to identify karst groundwater components [37,38]. We believe that δ 18 O is more appropriate for GMM, because rapid infiltration of heavy precipitation in our system can cause a change in groundwater EC in both directions, but only toward more positive values in the case of δ 18 O.
We used GMM to determine the probability densities of the δ 18 O groundwater series. GMM decomposes the total probability density into the sum of normal sub-densities (components). The contribution of each component to the total density was determined by a weighting factor according to relation (3). The analysis was based on the following assumptions: (a) a probability density composed of several sub-densities indicates the presence of water masses from different parts of the aquifer, e.g., from a porous matrix (base flow) and a conduit system (rapidly infiltrated water); (b) the standard deviation of the obtained Gaussian components comparable to the measurement error (0.1‰) is characteristic of a single water mass; and (c) the standard deviation larger than the measurement error is characteristic of component representing water mass with multiple contributions (e.g., a mixture of water from multiple precipitation events or a mixture of newly infiltrated water with a base flow).
The first part of the analysis covered the entire sampling period (April 2010-April 2012). The probability densities for W2 and DBC each have only one peak and the δ 18 O values are distributed over the normal distribution ( Figure 11, upper part). Corresponding standard deviations are larger than the measurement error of δ 18 O, indicating the presence of multiple water masses in the systems (Table 7).  To gain insight into aquifers functioning at the time of active recharge and discharge, we also considered the δ 18 O-probability density from the water year 2011 (October 2010-September 2011). In case of W2 for the period of water year 2011, we found the existence of two peaks (P1 and P2, Table 7, Figure 11). The mean value of sub-density corresponding to W2 s P1 for the water year 2011 is comparable to the δ 18 O values of groundwater samples collected during periods without significant precipitation and with low water level (Table 7, Figure 8), and the corresponding standard deviation is close to the value of the measurement error. These observations indicate that P1 represents the δ 18 O values of the base flow, i.e., water from the groundwater reservoir. The value of P2 is less negative than the value of P1 (Figure 11). The δ 18 O values of groundwater near the arithmetic mean of peak P2 occurred at the time of heavy rainfall and high flow (Table 7, Figure 8). Therefore, peak P2 can be attributed to the presence of rapidly infiltrated water. The weighting factors show the dominance of base flow (P1) relative to the component containing a contribution from rapidly infiltrated precipitation (P2) as it can be seen in Table 7.
This way of interpretation suggests that the δ 18 O values of groundwater change seasonally. During the colder part of the year there is heavy precipitation and vegetation activity is low. During this time, there is rapid infiltration of precipitation into the soil and activation of flow through wide karst channels, resulting in positive shifts in groundwater δ 18 O values. During the warmer part of the hydrologic year, high air temperatures and droughts occur. Groundwater δ 18 O values are more negative and fluctuate less than during the colder part of the year. This supports the conclusion that the base groundwater flow is present in MWs at this time, and that draining from the water reserves, which were primarily filled with winter precipitation and snowmelt, is taking place [25].
Using this way of reasoning, we conclude that baseline water was present in W2 during the dry periods of 2010 and 2011. The mean δ 18 O value for July-September 2010 (−8.64‰) is significantly lower than the mean groundwater value for May-October 2011 (−8.22‰) (t-test; t = −8.37; p < 0.001). The significantly higher values in 2011 can be attributed to the absence of significant snow events in the previous winter. Ultimately, this leads to the conclusion that the time of water storage in the subsurface is shorter than a year. This conclusion is consistent with the results obtained for this area [39]. Figure 8 shows that δ 18 O values on DBC are consistently higher than on W2, which is supported by statistical analysis (t-test; t = 3.29; p = 0.0006). One possible explanation is that DBC is recharged from lower altitudes than W2. However, DBC is very often influenced by the sea (chloride ion concentration: median = 22.8 mg/L, min = 8.4 mg/L; max = 420 mg/L). Because seawater is expected to have higher δ 18 O values than freshwater [16], we wanted to find out if seawater intrusion affected the isotopic content of DBC, and we used the MARS to test this hypothesis.

Multivariate Adaptive Regression Splines (MARS) for Estimating Isotopic Content of Dobrica (DBC) Spring
Statistical model MARS contained Cl (indicator of sea water presence), turbidity (indicator of newly infiltrated water), water temperature, and pH, as well as average air temperature and total precipitation between the two groundwater samplings. EC was not included in the analysis since it is strongly correlated with Cl (Table 6). Beside above mentioned parameters, the model also included the autoregressive component of first order of DBC's δ 18 O values (δ 18 O(−1)), since previous time-series analysis identified the AR(1) model as the best for DBC isotope content specification [25]. The obtained MARS equation (R 2 = 0.85) includes δ 18 O(−1), turbidity, precipitation, and Cl as significant determinants (Table 8). According to the MARS model, the autoregressive component has the greatest influence on the isotopic composition of DBC water, followed by turbidity (Table 8). Since the influence of the sea on the isotopic composition of DBC cannot be ignored, it is possible that the slightly higher values of the isotopic composition of DBC compared to W2 are the result of seawater intrusion into DBC, in addition to the different mean recharge altitudes of these two water resources. Figure 12 shows originally measured DBC δ 18 O values, results of the MARS model and residuals. The proportion of the variance in the dependent variable (DBC δ 18 O) that is predictable from the independent variables is 85%.

Conclusions
Karst aquifers are critical for drinking water supply in Mediterranean countries. They are extremely difficult to study and highly susceptible to pollution. Understanding how karst aquifers function is critical for proper management, sustainable use, and emergency response planning in the event of a major pollution event. Historically, hydrological discharge analyses and chemical-physical parameters have been used in the study of karst aquifers. This information, however, is not always available.
In the paper, we present an analysis of the activity of two karst groundwater resources, Martinšćica wells and Dobrica spring, which represent two adjacent karst aquifers in the Rijeka area (Croatia). Due to scarcity of traditional hydrological indicators, the statistical analysis of groundwater's stable isotopic composition received the most attention in the study.
The analysis revealed that precipitation from the cold part of the water year (October-April) is by far the most prevalent recharge to the systems, and the average residence time of water in the subsurface is less than a year. Stable isotopes were found to be a better indicator of the presence of newly filtered water in the system than the standard indicator of electrical conductivity. The Martinšćica wells system was shown to be a typical dualporosity system with dominant base flow. Using a nonparametric regression analysis, we demonstrated that the possibility of seawater intrusion into the spring affecting Dobrica's water isotope content cannot be ruled out. Nonetheless, the autoregressive component has the greatest influence on the isotopic composition of this spring.
We believe that the results presented in the paper are credible in demonstrating that a thorough statistical analysis of water's stable isotope content is not only an addition to traditional hydrological analysis, but is also appropriate for determining the hydrological behavior of karst systems when conventional parameters are unavailable or have limited availability.