Revealing hydrological relations of adjacent karst springs by partial correlation analysis

Hydrological and hydrogeological characteristics of karst aquifers are complex and significantly different from characteristics of hydrogeological environments such as fractured and granular aquifers (Bakalowicz 2005; Hartmann et al. 2012). Underground structures of fissures, pores, fractures and conduits of various size cause a high degree of heterogeneity, complex hydraulic conditions and spatial and temporal variability of the hydraulic parameters (Labat et al. 2001, 2016).

Preliminary studies (Padilla & Pulido-Bosch 1995; Larocque et al. 1998; Bouchaou et al. 2002; Panagopoulos & Lambrakis 2006; Mayaud et al. 2014) showed that utilizing of time-series analysis techniques (auto-correlation function (ACF) and cross-correlation function (CCF)) in karst hydrology gives a basic understanding of the functioning of the karstic system. The ACF represents the linear dependency of a time-series over a time period (Larocque et al. 1998). The CCF is used to determine the relationship between the input and the output time-series (Larocque et al. 1998; Panagopoulos & Lambrakis 2006). Hydrological time series are exhibited to diverse processes involved in the transfer of water in hydrological cycle. According to Jukić & Denić-Jukić (2011, 2015), effects of these processes can be displayed in ACFs and CCFs and consequently ambiguities with respect to the effects encoded in the correlation function exist.

The correlation coefficient quantifies the linear association between two time series (Sun et al. 2013). According to Liang (2014), correlation does not carry the needed directedness or asymmetry and therefore does not necessarily indicate causation. Correlation between two time series may be the result of another factor controlling both series (Jukić & Denić-Jukić 2015). The partial correlation coefficient measures the linear relationship between two time series whilst controlling for the effect of one or more (control) additional time series. While the first-order partial correlation coefficient measures the effect of a single control time series, the high-order partial correlation coefficient measures the total effect of several control time series. The number of control time series defines the order of partial correlation. Partial correlation coefficients have been applied in hydrology and water resources (Bellin & Rinaldo 1995; Vervier et al. 1999; Zhang et al. 2007; Burn 2008; Boucher et al. 2009; Vandenbohede et al. 2011; Fan et al. 2013), while partial cross-correlation function and matrix have been used mostly in other fields of science (Ehrentreich 1997; Stark 2006; Tkach et al. 2007; Enquist et al. 2008) until Jukić & Denić-Jukić (2011) applied partial correlation matrix for the first time to the field of hydrology. Jukić & Denić-Jukić (2015) proposed a new statistical method based on analysis of higher-order partial correlation functions for investigating relationships between rainfall and karst-spring discharge. The method showed improvement in simple correlation analysis because it overcomes the bound in the application of this analysis relating to ambiguities with respect to the impacts of space-time-variant processes encoded in time series.

This paper presents the application of simple correlation and first-order partial correlation analysis in order to study the behavior and relationship of two adjacent karst springs in a dry and wet hydrological year by using hourly data. The study site is the catchment of the karst springs Jadro and Žrnovnica in Croatia. This work focuses on the effects of temperature, relative humidity (RH), rainfall and discharge on the hydrological behavior and relations of two adjacent karst springs. According to Jukić & Denić-Jukić (2015), the groundwater recharge depends on the process of evapotranspiration and rainfall, so air temperature, relative humidity and rainfall control the spatial and temporal distribution of this process. This paper is the first that applied discharge as a control signal in a partial correlation analysis of karst springs and the first that applied this type of analysis on two adjacent karst springs.

The study site is the catchment of the karst springs Jadro and Žrnovnica that is geographically situated in Dalmatia (Croatia) (Figure 1). The catchment is located at the foot of Mosor (1,300 m a.s.l.) and Kozjak (600 m a.s.l.) mountains, mostly formed by carbonate rocks and partly formed by impermeable flysch.

The existing hydrological delineation is not accurately and reliably defined and described (Jukić & Denić-Jukić 2008, 2015; Kadić et al. 2015). Hypothetical boundaries cover approximately 450 km² (Figure 1). The north boundary represents a complete hydrogeological barrier located north of the Mućko Polje. The tracer tests of the ponor Jablan revealed a fast connection between this area and the springs Jadro and Žrnovnica. The catchment boundary on the west presents only a hypothetical delineation because its determination is based on less important hydrogeological indicators. The tracer tests of the ponor Postinje showed that the western part of the Mućko Polje is shared with catchments of the Ribnik Spring and the Mandrača Spring. The southern boundary of the catchment is the coastal barrier formed by the contact between permeable rocks of hinterland and impermeable rocks of the coastal region. The outlets of the Jadro and Žrnovnica Springs are situated in the lowest part of this contact area. The eastern boundary of the catchment is situated near the Cetina River that flows through a carbonate formation at elevations above 300 m a.s.l. (Jukić & Denić-Jukić 2008). The construction of hydroelectric power plants along the Cetina River has increased the discharge of the Jadro Spring and prevented the Žrnovnica Spring from drying up (Jukić & Denić-Jukić 2008). The tracer tests of the ponor Grabov Mlin that were performed in 1963 showed the direct connection between the Cetina River and the Jadro and Žrnovnica Springs (Bonacci 1987). The tracer appeared at both springs 23 days after injection into the ponor. All the results above indicate the existence of inter-catchment groundwater flows (Bonacci & Andrić 2015; Kadić et al. 2015).

The analyzed time period is the period of two hydrological years: (a) dry hydrological year from 20.10.2011 to 12.9.2012; and (b) wet hydrological year from 13.09.2012 to 09.09.2013. The available hourly data are:

• rainfall observed at meteorological stations Vučevica (V), Dugopolje (D) and Tunel Konjsko (TK),

• average air temperature (T) from meteorological stations Vučevica, Dugopolje and Tunel Konjsko,

• average relative humidity (RH) from meteorological stations Vučevica, Dugopolje and Tunel Konjsko,

• discharge of the Jadro Spring (J),

• discharge of the Žrnovnica Spring (Ž).

Basic statistical characteristics of time series data are presented in Table 1.

In karst catchment areas, the interpretation of hydrological, geological and hydrogeological research methods is limited due to the limited scope of such research. It has been shown that preliminary studies based on simple correlation analysis techniques can be performed to get a basic understanding of the functioning of the karst system (Padilla & Pulido-Bosch 1995; Larocque et al. 1998; Bouchaou et al. 2002; Panagopoulos & Lambrakis 2006).

The ACF quantifies the linear relationship of successive values within a signal, dependent on their time distance (Terzić et al. 2012). The CCF is used to analyse the linear relationship between an input and output signal (Zhou et al. 2014). Unfortunately, the CCF does not contain information indicating the cause of the linear relationship. Thus, the correlation of two signals may be the result of a third signal that controls or generates both analyzed signals (Jukić & Denić-Jukić 2011, 2015). This problem can be solved by using partial correlation analysis.

The partial correlation coefficient measures the linear relationship between two time series whilst controlling for the effect of one or more (control) additional time series.

Using a control signal that is uncorrelated with the input and the output signal (e.g. white noise signal), in Equation (7) the transition from the partial cross-correlation function to the cross-correlation function is evident: If r_xz → 0 and r_zy(k) → 0, then r_xy|z(k) → r_xy(k).

In the following analysis x signifies the input signal, y signifies the output signal and z signifies the control signal. The following abbreviations for functions are used:

• ACF: auto-correlation function;

• CCF: cross-correlation function;

• PACF: partial auto-correlation function;

• PCCF: partial cross-correlation function.

Interpretations of possible relations between the cross-correlation function CCF (x-y) and the partial cross-correlation function PCCF (x-y|z) are given in Figure 2.

Relations given in Figure 2 can be applied to the relations between auto-correlation and partial auto-correlation functions since the partial auto-correlation function is a special form of the partial cross-correlation function.

The ACFs of input rainfall signals and output discharge signals in the dry year are presented in Figure 3(a). The ACFs of hourly rainfall at all three meteorological stations diminish very rapidly. The rainfall signals are not auto-correlated after the time lag of 15 h at the Vučevica, 13 h at the Dugopolje and 9 h at the Tunel Konjsko. The memory effect, defined according to Mangin (1984) as the time lag where the ACF of spring discharge becomes less than 0.2, is 340 hours (about 14 days) for Jadro Spring and 324 hours (about 13 days) for Žrnovnica Spring. The truncation point is chosen to be 2,000 hours because ACFs and CCFs at lags larger than 2,000 hours contain only information about seasonal periodicity (Jukić & Denić-Jukić 2015).

Figure 3(b) presents the ACFs of input rainfall signals and output discharge signals in the wet year. The truncation point is chosen to be 2,000 hours as well as for the dry year. The auto-correlation functions ACF (V), ACF (D) and ACF (TK) decrease very quickly and reveal that the rainfall signals are not auto-correlated after the time lag of 5 h at the Vučevica, 7 h at the Dugopolje and 4 h at the Tunel Konjsko. The memory effect for Jadro Spring is 1,404 hours (about 58 days) and 1,260 hours (about 52 days) for Žrnovnica Spring.

Differences in memory effect for the dry and wet year (Figure 3(a) and 3(b)) are the consequence of the seasonal periodic component in discharge time series, which is mainly due to evapotranspiration. This component is more dominant in the wet year. Consequently, the memory effect for the Jadro Spring shows significant differences in relation to previous studies (Jukić & Denić-Jukić 2011, 2015) that applied auto-correlation analysis to longtime periods of daily data without taking into account the rainy and dry seasons. Correlation analysis was not applied to Žrnovnica Spring, therefore it is not possible to give a comparison.

In the dry hydrological year (Figure 4(a)) the cross-correlation functions CCF (J-V), CCF (J-D) and CCF (J-TK) have a shape similar to the cross-correlation functions CCF (Ž-V), CCF (Ž-D) and CCF (Ž-TK) (Figure 4(b)), although the discharge from Žrnovnica Spring shows better correlation with rainfall that is most evident in the peaks of the CCFs. The results of the cross-correlation response times are presented in Table 2. The response time is the time lag that corresponds to the maximum of the CCF (Delbart et al. 2016). According to Mangin (1984) the response time associated to the cross-correlation between rainfall and discharge corresponds to the mean response time of the karst aquifer to a rainfall event. In the case of the wet hydrological year (Figure 4(c) and 4(d)) the CCFs showed similar behavior for the Jadro and Žrnovnica Spring with a better correlation of the Jadro Spring to all meteorological stations excluding the first ‘peak’.

According to the results of cross-correlation response times presented in Table 2, both springs show rapid response (less than 24 hours) to a rainfall occurrence with almost equal values in the dry and the wet year. These results indicate the existence of an important quick-flow (Lo Russo et al. 2015) component for both springs. Statistically significant values at negative time lags are evident in the CCFs (Figure 4(a) and 4(d)). Since the rainfall-discharge relation is one-directional, these values confirm that the time series of rainfall and discharge contain the effects of other processes involved in the water transfer (Jukić & Denić-Jukić 2015).

The CCF between the Žrnovnica and Jadro Spring discharge (Figure 4(e)) results in high values and is mainly symmetrical in the dry and in the wet year, which indicates that the runoffs on both springs are generated by the same hydrological processes (Mayaud et al. 2014). Higher values on the positive side of the CCF in the dry and in the wet year imply the existence of a causal connection between the Žrnovnica Spring and Jadro Spring. For a detailed analysis of this connection partial correlation analysis was used.

Comparing the dry and wet year, both springs Jadro and Žrnovnica showed similar behavior. Both springs showed characteristics of an intermediate karstified aquifer (Padilla & Pulido-Bosch 1995).

The following first-order partial correlation analysis is focused to determine the influences of air temperature, RH, precipitation and discharges separately, on the hydrological behavior and relations of Jadro and Žrnovnica Springs in the dry and wet year. The time series of rainfall from the meteorological stations Vučevica (V) and Dugopolje (D) are treated as input signals. The time series of discharge of karst springs Jadro (J) and Žrnovnica (Ž) are treated as output signals. The control signals are time series of average air temperature (T), average RH, rainfall (TK) (that was chosen considering the location of the Tunel Konjsko meteorological station) and the discharges of Jadro and Žrnovnica Springs.

Figure 5(a) presents the comparison of the auto-correlation functionACF (J) and partial auto-correlation functions PACF (J|T), PACF (J|RH), PACF (J|TK) and PACF (J|Ž) in the dry year. It can be noted that the existing differences between ACF (J), PACF (J|RH) and PACF (J|TK) are insignificant, meaning that control signals RH and TK have no effect on the aquifer emptying. The effect of temperature (T) on the ACF (J) is somewhat more evident but still relatively low. Contrarily, PACF (J|Ž) differs significantly from ACF (J). Control signal Ž has partial explanation effect on the ACF (J) at legs less than 1,034 hours and suppression effect on the ACF (J) at lags 1,034–1,392 hours (Figure 5(a), Detail A). Similar results are obtained for the Žrnovnica Spring (Figure 5(b)) in the dry year where PACF (Ž|T) and PACF (Ž|J) showed divergence from ACF (Ž).

The effects of temperature (T) on the ACF (J) and ACF (Ž) in the wet year (Figure 5(c) and 5(d)) are more significant than in the dry year because the seasonal periodic component of discharge in the wet year is more significant as the differences in the amount of runoff during the summer and winter months are greater. Runoff in the dry year is more uniform over the year. This result shows that the inflows from the adjacent Cetina River catchment are more dominant during dry years, which is in accordance with previous studies (Stepinac 1983; Bonacci 1987) that showed that after the build of a series of accumulations on the Cetina River there was an increase in the minimum flow rates in the dry season for both springs, Jadro and Žrnovnica. It can be noted that control signals RH and TK have no effect on the ACF (J) and ACF (Ž) in the wet year (Figure 5(c) and 5(d)). Control signal J has a partial explanation effect on the ACF (Ž) for the first 340 hours (Figure 5(b), Detail A) and from 840 to 1,400 hours (Figure 5(b), Detail C), and suppression effect on the ACF (Ž) at lags 432–840 hours (Figure 5(b), Detail B). After 1,400 hours control signal J has no significant effect on the ACF (Ž).

Comparing the effects of PACF (J|Ž) and PACF (Ž|J), it is evident that the control signal J has more effect on ACF (Ž) than inversely in the dry and wet year (Figure 5). These results can be explained by the fact that the Žrnovnica Spring is partly recharged from the Jadro basin which is in accordance with the study from Bonacci & Andrić (2015). Bonacci & Andrić (2015) explained that probably after heavy rainfall events part of the Jadro basin drains into the Žrnovnica basin.

According to Jukić & Denić-Jukić (2015), the first break in the partial auto-correlation function denotes the end of the quick-flow and the second break denotes the end of the intermediate-flow and the actual system memory. There is no exact way to determine these breakpoints. The breakpoints are determined by noticing changes in the slope of the partial correlation function. The breakpoints in the partial auto-correlation functions in the dry year cannot be accurately determined (Figure 5(a) and 5(b)). In the wet year the breaks in the slopes of PACF (J|T) and PACF (Ž|T) are observable at the lag of 336 hours (14 days) (Figure 5(c) and 5(d)), which coincides with the results obtained by Jukić & Denić-Jukić (2015). It can be concluded that the quick-flow duration is 14 days for both springs. Unfortunately, other breaks in the slopes of partial auto-correlation functions PACF (J|T) and PACF (Ž|T) (Figure 5(c) and 5(d)) cannot be accurately determined. Partial correlation analysis was not applied to Žrnovnica Spring so far, therefore it is not possible to give a comparison.

Comparison of the cross-correlation function CCF (x-y) and partial cross-correlation functions PCCF (x-y|z) are presented in Figures 6–9. The interpretation of the results is related to the first 1,000 hours since the most significant effects of control signals occur during that period. Afterwards, the effects are mostly statistically insignificant.

The effects of the control signal T on the rainfall-discharge relations are observable but still relatively low (Figures 6 and 7), i.e. the partial cross-correlation functions PCCF (x-y|T) do not significantly differ from the cross-correlation functions CCF (x-y), regardless of the chosen input x and output signal y. Using relative humidity RH as the control signal resulted in partial cross-correlation functions PCCF (x-y|RH) that significantly differ from the cross-correlation functions CCF (x-y) (PCCF (x-y|RH) < CCF (x-y)) during almost the whole analyzed period of the dry year for all input and output signals (Figures 6 and 7). RH affects the process of evapotranspiration and thus the spatial and temporal distribution of the groundwater recharge process in the dry year (Jukić & Denić-Jukić 2011, 2015). Control signal TK has mostly a partial explanation effect on the rainfall-discharge correlation (Figures 6 and 7) in the first 1,000 hours.

The most significant suppression effect on the rainfall-discharge correlation in the dry year, which lasts for the first 48 hours, is obtained by using control signals J (Figures 6(a) and 7(a), Detail A) and Ž (Figures 6(b) and 7(b), Detail A). Afterwards, control signals J and Ž have almost no effect on the rainfall-discharge correlations (PCCF (x-J|Ž) ≈ CCF (x-J) and PCCF (x-Ž|J) ≈ CCF (x-Ž) for all input signals. The results of the partial cross-correlation analysis in the dry year are presented in Table 3.

Comparing PCCF (x-J|Ž) and PCCF (x-Ž|J) in the dry year (Figures 6 and 7), it is evident that the Žrnovnica Spring discharge has a greater effect on rainfall-discharge correlations than the Jadro Spring discharge, for all input signals.

In the wet year the effects of control signals T and RH are dominant almost during the whole analyzed period (Figures 8 and 9) for all input and output signals. Control signals T and RH remove the seasonal periodic component (Jukić & Denić-Jukić 2015) that comes out more in rainy periods, which means that the process of evapotranspiration significantly affects the groundwater recharge. Except for very short periods (especially the first 72 hours), control signal TK mostly suppresses the rainfall-discharge correlation (PCCF (x-y|TK) > CCF (x-y)) (Figures 8 and 9), i.e. it reduces the cross-correlation between signals x and y.

The suppression of the rainfall-discharge correlation in the wet year caused by control signals J (Figures 8(a) and 9(a), Detail A) and Ž (Figures 8(b) and 9(b), Detail A) lasts for the first 24 hours. After 24 hours control signals J and Ž have partial explanation effects on the rainfall-discharge correlations, i.e. PCCF (x-J|Ž) < CCF (x-J) and PCCF (x-Ž|J) < CCF (x-Ž) for all input signals. The results of the partial cross-correlation analysis in the wet year are presented in Table 4.

Comparing the partial cross-correlation functions PCCF (x-J|Ž) and PCCF (x-Ž|J) in the wet year (Figures 8 and 9), it is evident that Jadro and Žrnovnica Spring discharges have almost the same effect on the rainfall-discharge correlations, for all input signals.

According to Jukić & Denić-Jukić (2015), the first break in the partial cross-correlation function denotes the duration of the quick-flow and the second break denotes the intermediate-flow duration and the start of the pure base-flow. These results should be in accordance with results from the partial auto-correlation analysis. However, the breakpoints in the slopes of partial cross-correlation functions in Figures 6–9 cannot be accurately determined, therefore, no comparison for the Jadro Spring to the study from 2015 is possible.

The results of the partial correlation analysis (Figures 5–9) using discharge as a control signal show that a mechanism of groundwater exchange between the karst systems of springs Jadro and Žrnovnica exists, which is in accordance with the existing hypotheses about functioning of these two springs, i.e. these are springs with overlapping catchments that partly share the same aquifer.

The aim of this study was to investigate relations of two adjacent karst springs Jadro and Žrnovnica and their functioning during the dry and rainy seasons, applying simple correlation and first-order partial correlation analysis using hourly data.

The correlation analysis of karst springs Jadro and Žrnovnica showed similar behavior of both springs. Both springs showed characteristics of an intermediate karstified aquifer (Padilla & Pulido-Bosch 1995) with a much longer system memory in the wet hydrological year.

The results obtained by cross-correlation analysis showed that Žrnovnica Spring has better correlation to all meteorological stations in the dry year and Jadro Spring in the wet year. Rapid response of both springs (less than 24 hours) to a rainfall occurrence indicates the existence of an important quick-flow component for both springs.

The cross-correlation analysis between the Žrnovnica and Jadro Spring discharge resulted in extreme high values in the dry and in the wet hydrological year, which indicates the existence of a causal connection between the Žrnovnica Spring and Jadro Spring.

The results of partial auto-correlation analysis showed that RH and rainfall have no effect on the ACFs of Jadro and Žrnovnica Springs in the dry and in the wet year. The effects of temperature on the ACFs in the dry year are low but evident for both springs. In the wet year the effects of temperature are more significant in the case of both springs. Utilizing discharge as the control signal showed that although both discharges have great effect, the Jadro Spring discharge has more effect on the ACF of the Žrnovnica Spring discharge than inversely, in the dry and in the wet year. From the partial auto-correlation functions PACF (J|T) and PACF (Ž|T) the duration of the quick-flow for both springs is determined, in the wet year. The result for the Jadro Spring coincides with the result obtained by Jukić & Denić-Jukić (2015). Partial correlation analysis was not applied to Žrnovnica Spring so far, therefore it is not possible to give a comparison.

The partial cross-correlation analysis results showed that average RH affects the rainfall-discharge relation in the dry and in the wet hydrological year, unlike average temperature that has a significant effect only in the wet year. So it can be concluded that average RH affects the process of evapotranspiration in the dry year. In the wet year both average RH and temperature affect the process of evapotranspiration. Results of partial correlation analysis in the wet year showed that rainfall from the neighboring station can suppress the correlation between rainfall and karst discharge.

Applying discharge from adjacent spring as a control signal to partial cross-correlation analysis has proven that in the case of the springs Jadro and Žrnovnica, the discharge of each spring has an effect on the rainfall-discharge correlation and also suppresses the relation for the first 48 hours in the dry and 24 hours in the wet hydrological year.

Generally, the obtained results of correlation and partial correlation analysis show that these two adjacent karst springs share partly the same aquifer, which is in accordance with previous studies (Bonacci 1987; Jukić & Denić-Jukić 2008; Bonacci & Andrić 2015).

Partial correlation analysis represents an improvement in simple correlation analysis because it resolves ambiguities caused by influences of space-time-variant processes encoded in time series. For example, partial correlation analysis can be used for recognition of processes that affect the system response, determination of importance and contribution of each process, etc. Partial correlation analysis has two limitations in quantitative analyses: reliability and linearity. The reliability increases proportionally to the sample size. Since the partial correlation coefficient calculation is based on the simple correlation coefficient, the linear relationship is assumed. This limitation is not so important in qualitative analyses. The aim of qualitative analyses is to determine the behavior of the system and the reasons that cause such behavior. Exact quantitative values are not a priority in this type of analysis (Jukić & Denić-Jukić 2015).

The demonstrated approach represents an improvement in investigating adjacent karst springs in dry and wet periods, and can be applied on any adjacent karst springs. The proposed new approach could be used as an integral part of the hydrological and hydrogeological investigations of adjacent karst springs, considering that the necessary data can occur in a relatively simple and inexpensive manner.