An Application of Correlation and Spectral Analysis in Hydrological Study of Neighboring Karst Springs

Abstract

Various methods of time series analysis have been used in studies of karst hydrological systems. Among these methods, correlation and spectral analysis have had an important role. The correlation analysis most often has been based on determination of correlation coefficients and correlation functions. Partial correlation functions (PCF) are a mathematical tool of the correlation analysis which practical applicability in karst hydrology is insufficiently explored. In this study, the correlation and spectral analysis are applied on the catchment of Rumin Springs located in the Dinaric karst area between Croatia and Bosnia and Herzegovina. The available daily data are the rainfall, air temperature and relative humidity from three locations, as well as the discharge from two springs. The periods before and after the construction of HPP Orlovac in 1973 are analyzed. The basic hypothesis is that a difference between PCF obtained for two neighboring karst springs describe a difference in their functioning. The results of application show that PCF can resolve some ambiguities concerning the effects included in correlation functions and can provide the additional information that cannot be obtained by other methods of time series analysis. The obtained results are mostly in accordance with the present knowledge, and they support the existing hypotheses about the functioning of Rumin Springs.

1. Introduction

Karst is a terrain with distinctive landforms and drainage. A weathered zone of enhanced porosity named epikarst is situated at surface, while networks of fractures, fissures and conduits of various sizes and forms exist in underground, which all together represent a complex and very dynamic hydrological system [1,2,3,4]. The basic characteristic of karst terrains is a bimodal behavior. Concentrated turbulent flow takes place in conduits, while diffuse laminar flow prevails in fractures and fissures. These two hydrodynamic regimes can be also registered in karst spring discharge as quick flow and baseflow components. The network of conduits controls the quick flow component, whereas the network of fractures and fissures controls the baseflow component of discharge [5]. The existence of conduits makes karst terrains able to transport water rapidly on a regional scale. Consequently, karst terrains are extremely vulnerable to anthropogenic impacts [6], and they are sensitive to contamination from a variety of sources, so quantifying the proportion of the quick flow is often essential task. The consequence of internal structure of karst is the existence of fast fluctuations of groundwater level, usually with high amplitudes, producing quick alterations of velocity and direction of groundwater flow. It means that overflows from one to the other karst spring catchment occur depending on groundwater levels, so neighboring karst springs with overlapping catchments are common in many karst areas throughout the world [2,7,8]. The term “overlapping” illustrates the existence of partially separated systems, i.e., parts of catchment gravitating temporarily or permanently to several karst springs. The mechanism of groundwater exchange between neighboring karst springs is generally complex and time-variant [2,8,9,10,11].

Time series analysis comprises various methods for extraction of meaningful statistics and other characteristics of time series [12,13,14]. Among these methods, correlation and spectral analysis have had an important role in studies of karst hydrological systems. The correlation analysis most often has been based on determination of correlation coefficients, autocorrelation functions and cross-correlation functions. The spectral analysis has been based on determination of spectral density, cross-spectral density, phase, gain and coherence function [15,16,17]. These functions have found various applications in karst hydrology. They have been used for determination of basic characteristics of karst underground including investigations of groundwater hydrodynamics [18,19,20,21,22], processes on karst surface [23], transport properties [24,25], interactions between surface and underground flows [26], cavern conduit systems [27], connections in karst underground network [28], effects of storm behavior [29], residence and response times [30,31,32], transient inter-catchment flows [33], etc. In fact, correlation and spectral analysis have been used during last three decades mostly in preliminary hydrological analyses, where they can provide information supporting general conclusions about the functioning of karst hydrological systems [34,35,36,37,38,39].

Cross-correlation function is a measure of the strength of the linear relationship between two time series depending on a time lag between them. Autocorrelation function measures the strength of the linear relationship between successive values of a time series depending on a time lag between them. In addition to linearity, the correlation functions have also other drawbacks. Hydrological time series are affected by numerous space-time-variant processes occurring through water transfer inside hydrological cycle. The effects of these processes may be exhibited in the correlation functions. Consequently, in practical applications, it is very difficult to distinguish the contribution of each process, so ambiguities concerning the effects encoded in the correlation functions usually exist. Partial correlation functions sometimes can overcome these drawbacks. The partial cross-correlation function measures the remaining strength of the linear relationship between two time series after the removing the linear effect of other (control) time series. The partial autocorrelation function measures the remaining strength of the linear relationship between successive values of a time series after the removing the linear effect of other time series. They are a relatively new mathematical tool of correlation analysis in karst hydrology [40,41,42], which has been used for the analyses of the effects of meteorological and hydrological time series on the karst spring discharge [43,44,45], but their practical applicability generally is still insufficiently explored.

In this study, the correlation and spectral analysis are applied on the catchment of Rumin Springs. The Rumin Springs consist of two neighboring karst springs, named Rumin Veliki for the permanent one and Rumin Mali for the intermittent one. These two springs are situated in Croatia, but their catchment is shared between Croatia and Bosnia and Herzegovina. Partial cross-correlation functions are used to explain the existing relationships between the rainfall observed at three locations and the discharge from Rumin Mali and Rumin Veliki. Partial autocorrelation functions are used to estimate the effect of air temperature and relative humidity on the discharge from these two springs. This work has two aims: (1) advances in application of correlation analysis in hydrological studies of karst systems, with special reference to the hydrological study of neighboring karst springs by means of partial correlation functions, and (2) extension of present knowledge about the hydrological functioning of the Rumin Springs.

2. Materials and Methods

2.1. Study Area

The catchment of the springs Rumin Mali and Rumin Veliki is a part of the Cetina River catchment, which is situated in the Dinaric karst area (Figure 1). Hydrogeological information about these two springs and the catchment of Cetina River have been published in several papers [46,47,48,49,50,51,52]. A hypsometric and a geological map of study area, as well as a detailed location and cave maps of two springs can be found in [50]. Consequently, only a short summary of basic hydrological characteristics is presented hereafter.

The Cetina River is a permanent karst river of the Adriatic Sea basin located entirely within Croatian territory. The catchment covers an area of approximately 4000 km² [48]. It is a typical example of the river in deep karst region with topographically complex landscape, where the surface and subsurface catchment boundary do not coincide, so a precise delineation of catchment is very difficult. The Cetina River catchment mostly lies in the mountain region with peaks over 1500 m a.s.l., while the rest of it morphologically belongs to the areas of karst poljes (Kupreško Polje, Duvanjsko Polje, Glamočko Polje and Livanjsko Polje in Figure 1. Livanjsko Polje is considered as the world’s largest karst polje, which together with the area of Buško Blato covers 433 km². Climatologically, the north-eastern (Bosnian) part of the catchment has a continental climate with hot and dry summers and mild and wet winters, while the south-western (Croatian) part is under strong influence of the Adriatic Sea and the Mediterranean climate [53,54,55,56]. Consequently, the precipitation regime on the Cetina River catchment is complex and variable. The amount of average annual rainfall varies and hypsometrically differs from 1100 to 2700 mm with an average of about 1380 mm. The catchment area dominantly consists of permeable carbonate layers from Paleogene period. Quaternary formations of alluvial and marsh deposits are situated in larger flat areas like karst poljes and valleys [52]. The catchment is placed in the karst terrain formed of very thick layers of limestones and dolomites with intensive karstification and complex hydrogeological conditions. The whole area represents a typical karst terrain with specific karst forms like ponors and other speleological sites, so water from rainfall rather easily infiltrates. Consequently, the catchment is characterized by high, rapid infiltration with fast underground flows and by almost total absence of surface flow. The exceptions are the areas of the karst poljes, where intermittent and permanent surficial watercourses are formed on impermeable alluvial deposits and marls. Most of the groundwater flow occurs within limestones, in fissures, fractures and conduits of different and irregular shapes and varying sizes. The values of apparent groundwater flow velocities vary in a wide range between 0.5 and 24 cm/s [46], depending on the hydrological period and water table position. The slow movement of water is related with the dry season.

The Cetina River originates at about 380 m a.s.l. at the foot of Dinara Mountain. The total length of river is about 105 km, and it empties into the Adraiatic Sea near the town of Omiš. The river flow regime follows meteorological regime and displays maximum during period of autumn, a secondary one during spring, and a minimum in summer. Generally, the river is characterized by high seasonality in flow volumes due to the karst characteristic of terrain and the sudden inflow of either large amount of rain or melting snow from nearby mountains. The average annual flow is 105 m³/s, but flow during the summer ranges from 4 to 6 m³/s [48]. Morphologically, the Cetina River has a valley composed of narrow canyons incised in carbonate rocks and zones of lateral valley widening. The upper course of this river includes the Cetina springs zone, the Peruča Lake, a 1.5 km long narrow carbonate canyon, and the section through the Sinjsko Polje where the river has a morphology of a typical plain river. Several permanent and intermittent karst springs are located along the left side of the Peruča Lake and the riverbed through the Sinjsko Polje, which get their water from the carbonate rocks at contacts with the impermeable marl rocks. The most important are presented in Figure 1. The springs Rumin Mali and Rumin Veliki are among them. The distance between these two springs is 640 m. The spring Rumin Mali is located at higher altitude (326.8 m a.s.l.) in comparison with the Rumin Veliki (307.6 m a.s.l.). The results of tracing tests that were performed long time ago (between 1957 and 1961, [46]) show that all these springs are fed by water from the Livanjsko Polje through underground drainage systems of Dinara and Kamešnica Mountains. Concerning the Rumin Springs, the tracing tests of the ponor in Čaprazlije (697 m a.s.l.) revealed that the largest part of water from this area goes to the spring Dabar. Tracer appeared also on the springs Dragović, Peruča, Rumin Veliki, Rumin Mali and Cetina Springs. The tracing tests of Veliki Ponor (693 m a.s.l.) and Opaki Ponor (703 m a.s.l.) reveled that practically all water from this part of Livanjsko Polje goes to the spring Rumin Veliki. More than 94% of the injected tracer appeared at this spring. The remaining small quantities of tracer were observed on the springs Dabar, Peruča, Rumin Mali, Kosinac, Ruda Velika and Grab. The karst springs Ruda Velika, Grab and Kosinac drain dominantly the area of Buško Blato (696–703 m a.s.l.).

The Cetina River is a valuable ecological resource because of the existence of numerous rare ecosystems and landscapes. On the other hand, the water is used for drinking purpose, electricity production, industry, and irrigation, so the quantity and quality of water is influenced by human impacts. This river has very important role in the energy system of Croatia. Due to large hydro-energy potential, five hydroelectric power plants have been built in the catchment (Figure 1), so the natural hydrological regime of the Cetina River and nearby springs has been changed significantly. Concerning the Rumin Springs, the most important is the construction of hydropower plant (HPP) Orlovac including the accumulations Buško Blato and Lipa, which started with operation in 1973. The consequences of this construction were significant changes of minimum, maximum, and total discharges from two springs [49,50]. It is interesting, that the distribution of water towards two springs was changed in a different way [51]. The upper and smaller spring Rumin Mali was less affected by the construction than the lower and larger spring Rumin Veliki.

2.2. Available Data

Meteorological data were measured at two meteorological stations, MS Sinj and MS Livno, which are 32 km apart. MS Sinj is situated in Sinjsko Polje at altitude of 308 m a.s.l., whereas MS Livno is in Livanjsko Polje at altitude of 730 m a.s.l. The available data from these two stations cover only the period after construction of HPP Orlovac, from 1996 to 2018, where the meteorological data are complete and reliable. The analyzed daily time series from this period are the rainfall (P_S, P_L), air temperature (T_S, T_L) and relative humidity (RH_S, RH_L). Daily rainfall also was measured in the village of Bitelić Donji (P_B), at a rain gauge situated at the foot of Dinara Mountain (RG Bitelić Donji), but only during the period 1999–2008. The measurement of daily discharge from the spring Rumin Veliki (Q_RV) has been carried out since 1948 and from the spring Rumin Mali (Q_RM) since 1950. However, the analyzed time series are from the periods 1950–1972 and 1996–2018, where all data are complete. The basic statistical characteristics of the daily time series for three analyzed periods, 1950–1972, 1996–2018 and 1999–2008, are presented in

Table 1. Basic statistical characteristics of analyzed daily time series.

	1950–1972		1996–2018								1999–2008
	Rumin Mali	Rumin Veliki	Livno			Sinj			Rumin Mali	Rumin Veliki	Bitelić Donji
	QRM (m3/s)	QRV (m3/s)	PL (mm)	TL (°C)	RHL (%)	PS (mm)	TS (°C)	RHS (%)	QRM (m3/s)	QRV (m3/s)	PB (mm)
Min.	0	0.26	0	−15.6	15	0	−7.2	22	0	0.04	0
Mean	2.74	19.2	3.21	10.2	74	3.18	13.1	69	1.61	6.91	3.47
Max.	16.5	106	125.1	29.1	100	153.4	30.4	100	9.38	87.6	187.7
S. D.	3.09	20.1	8.06	12.1	8.03	8.53	19.7	7.63	2.19	12.0	9.40

.

The permanent karst spring Rumin Veliki had mean annual discharge of 19.2 m³/s before 1973. After the construction of HPP Orlovac, mean annual discharge is only 6.91 m³/s, or 36% of the previous. The intermittent karst spring Rumin Mali had mean annual discharge of 2.74 m³/s. After 1973, it amounts 1.61 m³/s, which is 59% of the previous. For both springs, the maximum measured discharge belongs to the period of natural conditions of flow in the Cetina River. It is 106 m³/s for the spring Rumin Veliki, and 16.5 m³/s for the spring Rumin Mali. The air temperature during period 1996–2018 ranged between −7.2 and 30.4 °C, at MS Sinj, and between −15.6 and 29.1 °C, at MS Livno. The mean air temperature for the entire period is 13.1 and 10.2 °C, respectively. The mean relative humidity at MS Sinj and MS Livno during the analyzed period was 69 and 74%, respectively. Although the altitudes of MS Sinj and MS Livno are different, these two stations have a similar rainfall regime. Specifically, the mean daily rainfall for MS Sinj and MS Livno for the period 1996–2018 are 3.18 and 3.21 mm, which corresponds to the mean annual rainfall of 1162 and 1172 mm, respectively. A more detailed hydrological analysis reveals that, during the same period, the annual rainfall ranged from 831 to 1686 mm at MS Sinj, and from 808 to 1796 mm at MS Livno. Correlation coefficient between daily rainfall at these two stations was 0.85. However, daily maximums observed in each year were not correlated. The rainfall observed at RG Bitelić Donji during the period 1999–2008 had a slightly different regime. The mean value and the maximum value are higher in Table 1, as well as the standard deviation, which show that the rainfall at this location is more abundant.

2.3. Methods

The hydrological functioning of the Rumin Springs before and after the construction of HPP Orlovac are compared and investigated by means of the correlation and spectral analysis. It involves the determination of autocorrelation functions, cross-correlation functions, and spectral density functions [12,13,14,15,16,17], as well as partial autocorrelation functions and partial cross-correlation functions [40,41,42]. A short theoretical explanation of these functions and their expressions are presented below.

Let x represent the input time series and y represent the resulting output time series from a hydrological system. The covariance function between series x and y, both of length n, is defined by

$$c_{xy}\left(k\right)=\frac{1}{n}{\displaystyle \sum }_{t=1}^{n-k}\left(x_t-{\mu }_x\right)\left(y_{t+k}-{\mu }_y\right); k=0,\pm 1,\pm 2,\dots ,m$$

(1)

where μ_x and μ_y are the means of x and y, respectively. The truncation point m determinates the domain of this function, i.e., the time interval in which the analysis is carried out. It is recommended that $m=n/3$ [14,15]. The cross-correlation function (CCF) between time series x and y is

$$r_{xy}\left(k\right)=\{\begin{array}{c}\frac{c_{xy}\left(k\right)}{{\sigma }_x{\sigma }_y};k\ge 0\\ \frac{c_{yx}\left(-k\right)}{{\sigma }_x{\sigma }_y};k<0\end{array}$$

(2)

where σ_x and σ_y are the standard deviations of x and y, respectively. It should be noted that the correlation coefficient (CC) between time series x and y can be obtained as $r_{xy}=r_{xy}\left(0\right)$. For completely random time series, 95% confidence limits are approximately $\pm 2/\sqrt{n}$ [14]. The autocorrelation function (ACF) of these two series can be obtained from Equation (2) as

$$r_{xx}\left(k\right)=\frac{c_{xx}\left(k\right)}{{\sigma }_x^2}, r_{yy}\left(k\right)=\frac{c_{yy}\left(k\right)}{{\sigma }_y^2}$$

(3)

where $k\ge 0$. Spectral analysis is complementary to the correlation analysis. In fact, spectral density function (SDF) represents a Fourier transform of ACF. For $r_{xx}\left(k\right)$ and $r_{yy}\left(k\right)$, they are

$$S_x\left(f\right)=2\left[1+2{\displaystyle \sum }_{k=1}^{m}D\left(k\right)r_{xx}\left(k\right)\cos \left(2\pi fk\right)\right], S_y\left(f\right)=2\left[1+2{\displaystyle \sum }_{k=1}^{m}D\left(k\right)r_{yy}\left(k\right)\cos \left(2\pi fk\right)\right]$$

(4)

for frequencies $f=l/2m$, $l=0,1,2,\dots ,m$, where

$$D\left(k\right)=\frac{1+cos\frac{\pi k}{m}}{2}$$

(5)

represents the lag window proposed by Tukey [14,15], which ensures that the estimated values are not biased.

Let z represent a time series that controls the process of transformation of input series x to output series y. The linear effect of control series z can be removed from CCF, $r_{xy}\left(k\right)$, by calculating the partial cross-correlation function (PCCF):

$$r_{xy|z}\left(k\right)=\frac{r_{xy}\left(k\right)-r_{xz}{r}_{zy}\left(k\right)}{\sqrt{\left(1-r_{xz}^2\left(k\right)\right)}\sqrt{\left(1-r_{zy}^2\left(k\right)\right)}}$$

(6)

where $r_{xz}$ is the CC between series x and z, whereas $r_{zy}\left(k\right)$ is the CCF between series z and y [41,42]. The confidence limits are approximately the same as for CCF; $\pm 2/\sqrt{n}$. The linear effect of control series z can also be removed from ACF, $r_{xx}\left(k\right)$ or $r_{yy}\left(k\right)$. The result is a partial autocorrelation function (PACF). For example, PACF for $r_{yy}\left(k\right)$ is

$$r_{yy|z}\left(k\right)=\frac{r_{yy}\left(k\right)-r_{yz}{r}_{zy}\left(k\right)}{\sqrt{\left(1-r_{yz}^2\left(k\right)\right)}\sqrt{\left(1-r_{zy}^2\left(k\right)\right)}}$$

(7)

The effect of control time series z is investigated by comparing CCF and PCCF as well as ACF and PACF. Four types of effects of control time series z can be distinguished [42]:

• no effect—input time series x is an antecedent cause of output time series y and there is no relationship between series z and y,
• explanation effect (control effect)—control series z is an antecedent cause of output series y and there is no relationship between series x and y,
• partial explanation effect—input series x is an antecedent cause of output series y, but control series z is also either antecedent or intervening cause,
• suppression effect—control series z has a positive effect through one path and a negative effect through another path, or when control series z is weakly correlated with one of the original series.

Generally, the values of PCCF may be the same, lower, or higher than CCF. For example, Figure 2 presents schematically the effects for a hydrological system with two input time series of correlated rainfall, x and z, and an output time series of discharge y. The following general rule for the notation of functions is used: $\mathrm{NAME}(x,y|z)$, where $\mathrm{NAME}$ is a function abbreviation (CCF, ACF, PCCF or PACF), x is an input time series, y is an output time series, and z is a control time series.

3. Results

The autocorrelation functions of discharge Q_RM and Q_RV for periods 1950–1972 and 1996–2018 are presented in Figure 3. Generally, ACF quantifies the linear dependency of successive values over a time and outlines a system memory. In karst hydrology, ACF of spring discharge provides the information about the storage capacity of karst system. The storage capacity depends on the degree of karstification, but also on the surface ponds, soil cover and epikarst properties. To compare functions of various aquifers, the so-called memory effect is defined as the time lag where ACF becomes less than 0.2 [15]. An undeveloped karst system, with a large storage, has a high memory effect where ACF shows a slightly decreasing slope. On the contrary, a developed karst network, without important storage, corresponds to a low memory effect where ACF has a much steeper slope. In Figure 3, it can be noted that ACF(Q_RM) and ACF(Q_RV) for the first period 1950–1972 have very similar forms and the memory effect of approximately 63 days. Concerning the second period 1996–2018, the form of ACF(Q_RM) is only slightly changed, as well as the memory effect that is reduced for only 1 day, and it amounts 62 days. It shows that, despite an evident change in hydrological regime after construction of HPP Orlovac in 1973, the storage capacity of the spring Rumin Mali apparently remains the same. It is not the case with the spring Rumin Veliki. The values of ACF(Q_RV) obtained for the period 1996–2018 are significantly lower, and the memory effect is reduced to approximately 40 days. Concurrently, a bimodal behavior of this spring has become more evident. It should be emphasized that the time series of karst spring discharge contain a seasonal oscillation, so ACF(Q_RM) and ACF(Q_RV) in Figure 3 exhibit an evident periodicity with a period of approximately 365 days. The seasonal oscillation of discharge in Dinaric karst is dominantly the consequence of a similar oscillation in evapotranspiration rate during a year, which is manifested in groundwater recharge.

The spectral density functions of discharge Q_RM and Q_RV for the period 1950–1972 and 1996–2018 are presented in Figure 4a–d. SDF shows how the variance of analyzed time series is distributed over the range of frequencies. The peaks in SDF at various frequencies lead to the identification of periodicity in series. In system analyses, SDF of output time series also contains information about the filtering effect of analyzed system. It implies that there is a possibility to identify a change in the filtering effect during a time by comparing SDF of output time series obtained for two periods. In this study, SDF is used to analyze the change in short period and long period oscillations of discharge from two analyzed springs. In this context, the term long period implies the oscillations with the periods larger than 10 days, which corresponds to the frequency range below 0.1 1/day (Figure 4a,b). Truncation point is $m=2556 \mathrm{days}$ (approximately 7 years). It can be noted that the seasonal oscillation of karst spring discharge is manifested in both SDF(Q_RM) in Figure 4a, and both SDF(Q_RV) in Figure 4b, as a periodicity at frequency 0.00274 1/day corresponding to the period of 365 days. The comparing SDF(Q_RM) for two periods in Figure 4a reveals that the variance associated to the seasonal periodicity in discharge Q_RM is only slightly attenuated for the second period. The attenuation of the seasonal periodicity is more evident for discharge Q_RV (Figure 4b), which also indicates that the storage capacity of the aquifer of Rumin Veliki is significantly lower after the human intervention in 1973. Concerning the short period oscillations of discharge with the periods less than 10 days, which corresponds to the frequency range above 0.1 1/day, the obtained results are presented in Figure 4c,d. The truncation point is $m=365 \mathrm{days}$. This range mostly describes oscillations in peak discharge. It can be noted that the values of SDF(Q_RM) for the period 1950–1972 and 1996–2018 are similar (Figure 4c). Two functions have only small differences in shape, which means that the distribution of variance for the discharge from Rumin Mali has not been significantly changed at high frequencies. Concerning the discharge from Rumin Veliki (Figure 4d), the distribution of variance for two periods is evidently different. The values of SDF(Q_RV) for the period 1996–2018 are higher, which shows that the filtering effect of aquifer is reduced. In fact, a redistribution of variance from the low to high frequency range occurred after 1973, so the relative importance of the oscillations in peak discharge increased.

The spectral density functions for time series of rainfall P_L and P_S obtained for the period 1996–2018 are compared in Figure 5a. It is evident that MS Livno and MS Sinj have a very similar regime of daily rainfall. They both have similar periodicities and a seasonal component revealing that time series P_L and P_S cannot be considered as a white noise. A similar result is obtained also for RG Bitelić Donji for the period 1999–2008 (Figure 5b). The seasonal periodicity in rainfall is a climatological characteristic of study area.

The cross-correlation functions between rainfall (P_L, P_S) and discharge (Q_RM, Q_RV) for the period 1996–2018 are presented in Figure 6, as well as a confidence interval. CCF is an indicator of linear dependency between two time series as a function of time lag k. If CCF shows statistically significant values and it is not symmetrical, a causal relationship between two series exists. In a system analysis, if the input time series is a white noise, CCF represents the impulse response function. If the input time series contains a periodicity, CCF may exhibit an oscillation depending on the filtering effect of the system. Considering a karst hydrological system, CCF between the input time series of rainfall and the output time series of discharge provides the information about the system response including significance, duration, and time delay. The duration of the system response represents the range of positive lags from the origin where CCF has statistically significant values. The time delay (time to peak or response time) gives an estimate of the pressure pulse transfer times through the karst system. It also indicates the degree of karstification. Consequently, several details require attention in Figure 6:

• The seasonal periodicity observed in rainfall (Figure 5a) is manifested in all functions, especially at negative lags where statistically significant values exist. This periodicity is more evident for the discharge from Rumin Mali (CCF(P_S,Q_RM) and CCF(P_L,Q_RM) in Figure 6) than for the Rumin Veliki (CCF(P_S,Q_RV) and CCF(P_L,Q_RV) in Figure 6). This result is in accordance with the results of spectral analyses, where the spectral density function for the Rumin Mali in Figure 4a obtained for the period 1996–2018 has larger values at frequency 0.00274 1/day than the spectral density function for the Rumin Veliki in Figure 4b for the same period.
• Differences between CCF(P_S,Q_RM) and CCF(P_L,Q_RM) and between CCF(P_S,Q_RV) and CCF(P_L,Q_RV) are small, which shows that the rainfall from MS Sinj and MS Livno produce similar cross-correlation functions. The same result is obtained also for RG Bitelić Donji for the period 1999–2008, which means that the entire catchment probably has a similar variation of daily rainfall, considering a long period.
• The time delay is the same for all functions. It amounts 1 day, so the pressure pulse is quickly transferred towards both springs after rainfall.
• CCF(P_S,Q_RV) and CCF(P_L,Q_RV) have larger values at lags below 10 days than CCF(P_S,Q_RM) and CCF(P_L,Q_RM), which means that the discharge from Rumin Veliki is better correlated with rainfall during the period of quick flow. On the other hand, the discharge from Rumin Mali is better correlated with rainfall during the period of baseflow. The spring Rumin Mali also has evidently longer duration of response. Specifically, CCF(P_S,Q_RM) and CCF(P_L,Q_RM) become statistically insignificant at legs above approximately 135 days, whereas CCF(P_S,Q_RV) and CCF(P_L,Q_RV) become insignificant at legs above 92 days. These results show that a difference in the degree of karstification between the aquifers of two springs exist.

The partial cross-correlation functions between the input time series of rainfall P_L and the output time series of discharge Q_RM and Q_RV, where the control time series are rainfall P_S and P_B, are compared in Figure 7a,b. The main aim of this analyses is to characterize groundwater connections between the Livanjsko polje and the Rumin Springs. It should be emphasized that RG Bitelić Donji is situated at the foot of Dinara Mountain, immediately next to the catchment boundary of Rumin Mali (Figure 1). On the other hand, MS Sinj is in Sinjsko Polje, almost 10 km away from the catchment of Rumin Springs. It means that the control time series P_S theoretically has no effect (Figure 2). In practice it is not the case, because of the spurious correlation that always exists between rainfall observed at neighboring locations. The rainfall time series P_S evidently has a similarity with rainfall on the Dinara Mountain, so the effect of this similarity is removed from cross-correlation functions in Figure 6 by using P_S as the control series. The obtained results show that the seasonal periodicity is not present in PCCF(P_L,Q_RV|P_S), PCCF(P_L,Q_RM|P_S), PCCF(P_L,Q_RV|P_B) and PCCF(P_L,Q_RM|P_B). It is completely removed by using the control time series P_S and P_B, which reveals that the origin of this periodicity in cross-correlation functions is not the rainfall P_L. This result indicates that a relationship between the rainfall from Livanjsko Polje and the baseflow components of Rumin Springs does not exist. At lags corresponding to the first 10 days of spring response, PCCF(P_L,Q_RV|P_S) and PCCF(P_L,Q_RV|P_B), have evidently larger values than PCCF(P_L,Q_RM|P_S) and PCCF(P_L,Q_RM|P_B). In addition, control time series P_B explains almost completely the existing correlation between P_L and Q_RM, so PCCF(P_L,Q_RM|P_B) in Figure 7b fluctuates at upper boundary of confidence interval. It shows that the relationship between time series P_L and Q_RV is relatively strong, whereas the relationship between P_L and Q_RM hardly exists. This result indicates that the part of catchment in Livanjsko Polje has an important role in generation of the quick flow component of discharge from the spring Rumin Veliki. It is not the case with the spring Rumin Mali. Concerning lags larger than 10 days, the functions obtained for two springs have similar behavior. Statistically significant effects are registered approximately during first 60 days. Figure 7c,d compare the partial cross-correlation functions between the input time series of rainfall P_S and P_B and the output time series of discharge Q_RM and Q_RV where the control time series is rainfall P_L. The main aim of this analysis is to estimate the relative importance of the part of catchment located on the west side of Dinara Mountain. It can be noted that the partial effects in Figure 7c,d are more significant than the effects in Figure 7c,d, especially during the first 10 days corresponding to the quick flow. The strongest relationship at these lags is obtained for the rainfall from MS Bitelić (Figure 7d), which is in accordance with the present knowledge that the origin of quick flow components of both springs is the catchment on Dinara Mountain. The shapes of PCCF(P_B,Q_RV|P_L) and PCCF(P_B,Q_RM|P_L) are similar, which means that the responses of two springs are rather synchronized during the first 10 days after intensive rainfall. Concerning the partial effects located on other lags, results show that PCCF(P_S,Q_RV|P_L) and PCCF(P_S,Q_RM|P_L) in Figure 7c as well as PCCF(P_B,Q_RV|P_L) and PCCF(P_B,Q_RM|P_L) in Figure 7d consist the seasonal periodicity that was noted in CCF(P_S,Q_RV) and CCF(P_S,Q_RM) (Figure 6). The control time series P_L practically has no effect on this periodicity, which confirms that the catchment located on Dinara Mountain is mainly responsible for the generation of the baseflow components of both springs.

The effects of meteorological processes in catchment on the discharge from each spring are studied by means of the partial autocorrelation functions. The subjects of interest are the effects of rainfall P_L, air temperature T_L and relative humidity RH_L observed at MS Livno, on the form of ACF(Q_RM) and ACF(Q_RV). The obtained results are presented in Figure 8a,b. It can be noted that the control time series of rainfall P_L practically has no effects on the form of autocorrelation functions. The effects of control time series of relative humidity RH_L are observable, but they also do not play any important role. The largest effect is produced by the control time series of air temperature T_L. It should be noted that very similar results are obtained by using the meteorological data from MS Sinj, so they are not presented here. The resulting function for the spring Rumin Veliki (PACF(Q_RV|T_L) in Figure 8a) can be divided in two segments. First segment comprises lags below 60 days, where PACF(Q_RV|T_L) has statistically significant values. It shows that the actual memory effect of this spring is 60 days, which is 20 days longer than the estimation obtained from ACF(Q_RV) using classical method [15]. The second segment includes lags above 60 days where PACF(Q_RV|T_L) is rather flat, and it variates around the interval of confidence. This result shows that the time series T_L removes almost completely the effect of seasonal periodicity from ACF(Q_RV), which means that the air temperature has very important role in generation of the baseflow component of Rumin Veliki. Concerning the spring Rumin Mali, two segments also can be recognized in PACF(Q_RM|T_L) in Figure 8b. First segment covers lags below 63 days where PACF(Q_RM|T_L) has statistically significant values. The second segment includes lags above 63 days where, differently from the spring Rumin Veliki, time series T_L cannot remove periodicity. The existence of a seasonal periodicity in PACF(Q_RM|T_L) and the absence in PACF(Q_RV|T_L), as well as the evident differences between PACF(Q_RM|T_L) and PACF(Q_RV|T_L) at lags below 60 days, indicate that the springs Rumin Mali and Rumin Veliki have significantly different hydrological functioning.

4. Discussion

Generally, it is very difficult to explain reliably hydrological connections between neighboring karst springs and their catchments without extensive hydrological and hydrogeological investigations. The springs Rumin Mali and Rumin Veliki are two neighboring karst springs that are insufficiently explored. A precise catchment boundary delineation does not exist, and probably it is not possible due to catchment overlapping. Groundwater level monitoring has not been performed. There are only two meteorological stations in this area with a permanent monitoring, located in the Livanjsko and Sinjsko Polje. The significant drawback is a lack of information about the rainfall regime on the Dinara Mountain. According to the present knowledge [49], it seems that the springs Rumin Mali and Rumin Veliki are not hydrologically connected, even though their areal distance is only 640 m. Even if they have no significant groundwater exchange between catchments, similarity in hydrological regime of these two karst springs is evident, which is most often in karst the consequence of similar geological and climatological characteristics of catchments. The largest part of the catchment of Rumin Springs is in the Dinara Mountain. Both karst springs are recharged by the point and diffuse recharge from this mountain. The central part of this massif is an area of deep karst, where the depth of vadose zone amounts almost 1 km. However, due to highly karstified terrain, precipitation easily penetrates and feeds the deep karst aquifer. Because of deep karstification of carbonate bedrock and high amount of annual precipitation, the Dinara Mountain presents an important groundwater reservoir [57,58]. The present hydrological role of the part of catchment located in Livanjsko Polje is not completely clear. The results of tracing tests [46] show that probably, in its natural state before 1973, the spring Rumin Mali was supplied from the south-east part of Livanjsko Polje in very small quantities, in comparison with the quantities of water that fed the spring Rumin Veliki. Due to the construction of HPP Orlovac and accumulation Buško Blato, the recharge from the south-east part of Livanjsko Polje is strongly decreased, and now it probably exists only as occasional overflow, which is activated after intensive precipitation and flooding in some parts of Livanjsko Polje [50]. The relative importance of this overflow for the spring Rumin Veliki is unknown. However, the spring Rumin Mali is probably supplied in larger quantities than Rumin Veliki from the north-west part of Livanjsko Polje [49], which is not so much affected by the construction of HPP Orlovac. Consequently, hydrological analyses show that the performed interventions in the karst system of two springs had a lower influence on the hydrological regime of Rumin Mali than it had on the regime of Rumin Veliki [49,50,51].

The results obtained by time series analysis are in accordance with present knowledge about the consequences of construction of HPP Orlovac, and about the hydrological functioning of Rumin Springs after 1973. The comparations of the autocorrelation functions and the spectral density functions of the discharge from Rumin Springs for the period 1950–1972 and 1996–2018 reveal that the storage capacity of the aquifer of Rumin Mali is almost the same. At the same time, the storage capacity of the aquifer of Rumin Veliki is evidently smaller, and general characteristics of karst network apparently changed, which confirms that the catchment area of this spring was significantly reduced. The south-east part of the Livanjsko Polje was practically excluded from the catchment in 1973, so the relative importance of the highly karstified part of catchment located on the Dinara Mountain was increased, especially for the spring Rumin Veliki. However, it should be noted that the autocorrelation functions and spectral density functions of discharge also show that the storage capacities of the aquifers of Rumin Mali and Rumin Veliki were very similar before 1973.

The spectral density functions of rainfall from the meteorological stations MS Livno and MS Sinj obtained for the period 1996–2018 show that a seasonal periodicity in rainfall is present, which explains the existence of a similar periodicity in the cross-correlation functions between rainfall and discharge. Concerning the quick flow component, the results show that the discharge from Rumin Veliki is better correlated with rainfall than the discharge from Rumin Mali, which indicates a more developed karst network in this aquifer. However, the time delay is the same, 1 day for both springs. Concerning the baseflow component, the discharge from Rumin Mali is better correlated with rainfall, and it has also a longer duration of response, which also demonstrate that the aquifer of Rumin Mali has higher storage capacity than the aquifer of Rumin Veliki. The comparation of spectral density functions obtained for MS Livno and MS Sinj confirms that the variations of daily rainfall in Livanjsko and Sinjsko Polje are very similar. This similarity produces a spurious correlation between the daily rainfall observed at these two locations. The consequence of the spurious correlation is an extreme similarity between the cross-correlation functions obtained by using rainfall from MS Livno and MS Sinj as input series, which complicates interpretation of results.

The partial cross-correlation functions have been used to resolve the similarity between the cross-correlation functions, and to identify the difference in the responses of two springs on the rainfall events observed at MS Livno and MS Sinj. In this analysis, the daily rainfall from RG Bitelić Donji has been also used. This rain gauge is situated at the foot of Dinara Mountain immediately next to the catchment of Rumin Mali, and the rainfall from this rain gauge represents the rainfall regime on Dinara Mountain better than the rainfall from MS Sinj. It should be emphasized that the partial correlation analysis in this study is based on the comparation of the differences between the forms of partial cross-correlation functions obtained for two springs. The basic hypothesis is that a difference between the forms of partial cross-correlation functions obtained for two springs describe a difference in the functioning of these two springs. Because of the existence of significant spurious correlations between rainfall observed at neighboring locations, the rainfall data collected at locations around catchment boundaries could be used. Concerning the quick flow component of discharge, the obtained partial cross-correlation functions show that the catchment located in the Livanjsko Polje still has a role in generation of the discharge from the Rumin Veliki. A groundwater connection between the Livanjsko Polje and the discharge from the Rumin Mali cannot be confirmed reliably, because of the existence of a weak partial correlation in Figure 7b. Concerning the baseflow component, any significant connection between the Livanjsko Polje and the Rumin Springs has not been identified in this analysis. Generally, the results show that the groundwater reservoir of the Dinara Mountain has a crucial role in generation of the quick flow component and the baseflow component of both springs after 1973. The existing similarities between the partial cross-correlation functions obtained for two springs for the period of quick flow in Figure 7d show that pressure pulse transfer times are similar, which opens a possibility that two catchments are hydrologically connected during high water levels.

The partial autocorrelation functions of discharge have been used to estimate the effect of the meteorological processes in catchment on the discharge from each karst spring. The obtained results confirmed that rainfall and relative humidity have no significant effects on the form of the autocorrelation function of karst spring discharge [40,41]. However, the effect of air temperature is very important. Using the air temperature as the control time series, the seasonal periodicity has been reduced (case of the Rumin Mali) or almost completely removed (case of the Rumin Veliki) from the autocorrelation function of discharge, which indicates that the process of evapotranspiration has very important role in the generation of discharge, especially for the spring Rumin Veliki. The seasonal variation in discharge of Rumin Veliki mainly coincides with the seasonal variation of evapotranspiration. The evident differences in partial autocorrelation functions obtained for the springs Rumin Mali and Rumin Veliki demonstrate that these two springs have significantly different hydrological functioning, which indicates that the springs Rumin Mali and Rumin Veliki are supplied dominantly from two different catchments, at least during baseflow.

5. Conclusions

The correlation and spectral analyses are classical methods of time series analysis that have been often used in karst hydrology during last three decades. The correlation analysis usually has been based on the determination of cross-correlation coefficients, autocorrelation functions and cross-correlation functions. The partial correlation functions have been used rarely. However, the results presented in this work show that the partial autocorrelation function and partial cross-correlation function can be useful in the study of karst hydrological system functioning, especially in the situation of scarce data and information. By applying these functions, it is possible to resolve some ambiguities in the correlation functions and to provide the additional information that cannot be obtained by using other methods of time series analysis. Particularly, the partial cross-correlation functions have been used to resolve the similarity between the cross-correlation functions, and to identify the difference in the responses of two springs to rainfall in catchment. The partial autocorrelation functions have been used to resolve the seasonal periodicity of the autocorrelation functions, and to estimate the role of the process of evapotranspiration in the generation of the discharge from two spring.

The obtained results of the application of correlation and spectral analysis on the Rumin Springs are in accordance with the present knowledge, and they mostly support the existing hypotheses about the functioning of these springs. Particularly, it has been confirmed that the springs Rumin Mali and Rumin Veliki drain dominantly different parts of the catchment of Cetina River, which is manifested as a different hydrological functioning of two springs. The part that belongs to the spring Rumin Veliki has been much more affected by the construction of HPP Orlovac than the part that belongs to the spring Rumin Mali, which is exhibited as a significant reduction in the storage capacity of the aquifer of Rumin Veliki, whereas the storage capacity of the aquifer of Rumin Mali practically remains the same. It has been confirmed that the existing groundwater connection between the Livanjsko Polje and the spring Rumin Veliki affects the present hydrological regime of this spring, which can be explained with the occasional overflow activated after intensive precipitation and flooding in some parts of Livanjsko Polje.