A Research on Discharge and Temperature Regime of Substantially Altered Karst River

1. Introduction

The basic characteristic of karst is the existence of underground voids in limestone/dolomite rock mass. Karst underground voids are subsurface cavities formed primarily by the dissolution of soluble rocks (limestone, dolomite, gypsum, or salt). Following standard karst geomorphology and hydrogeology literature [1,2,3,4,5], the classification of karst underground voids can be approached by origin, geometry and morphology, hydrogeological function, and scale. From the standpoint of groundwater circulation, most important thing is scale. The scale implies a size-based classification where four types of voids are distinguished: (1) microvoids (dissolved pores and microfractures, < 1 mm), (2) mesovoids (enlarged fractures and small fissures, 1 mm – 10 cm), (3) macrovoids (conduits, channels and small caves, 10 cm – 10 m), and (4) megavoids (large, interconnected cave networks > 10 m). The size is controlled by lithology (solubility and purity of rock mass), structural geology (presence of faults, bedding and joints), hydrology (recharge type, base level and flow regime), chemical factors (concentration of CO₂ and acids in water), and climate. The size, geometry, and connectivity of voids directly control how groundwater flows, stores, and responds within the karst aquifer [6,7], so it defines the shape of karst spring hydrographs [8]. Diffuse, conduit and mixed flow through the karst aquifer can be distinguished [9,10,11,12]. The diffuse flow occurs in small voids such as microvoids so dominates in poorly developed or low-karstification zones (e.g., dolomitic terrains). It is slow, laminar flow following Darcy’s law. Recharge is diffuse, spread over a wide area. On a karst spring hydrograph, diffuse flow has several characteristic manifestations related to timing, shape, and behavior of baseflow recession. Diffuse flow dominance is shown by smooth rising limbs, moderate and broad peaks, long and gentle recessions, stable baseflow, and low flashiness. During recession period, temperature and chemistry of spring water are stable implying long residence time. The conduit flow is turbulent flow dominated by gravity and pressure gradients that occur in macrovoids. In addition to diffuse recharge through microvoids, conduit flow is formed also by concentrated recharge through vertical shafts. These voids drain the system quickly after recharge events, producing flashy hydrographs at karst springs. Hydraulic response to rainfall is fast, with short residence time so temperature and chemistry of spring water are unstable. The mixed flow exists in mesovoids with highly variable hydraulic conductivity depending on fracture and fissure density, aperture, and connectivity, where water flows through a combination of laminar and transitional regimes.

Consequently, water temperature is a valuable natural tracer of groundwater circulation and residence time in karst systems [13,14,15]. Stable temperature of karst spring water denotes long residence of infiltrated rainwater in underground karst system. Sudden temperature changes produced by influx of surface rainwater or stream denotes short residence time. Karst spring water temperature helps infer recharge events, flow path lengths, velocity, and connectivity between groundwater and spring [16,17,18,19,20]. Namely, infiltrating water often has different temperatures than aquifer water, leading to drops or rises in spring temperature during and after rainfall events. Along the flow path within a karst conduit, the thermal signal is changed by mixing, heat exchange, and storage effects, so temperature patterns provide spatial information on flow path length, conduit geometry, and exchange zones. Longer flow paths, contact with rock, and hyporheic zones reduce amplitude of temperature changes and introduce lags [21,22,23,24,25,26]. If water moves rapidly in large conduits, channels and caves temperature fluctuations will be large and rapid because of little contact between rock and water. If water flows slower through fractures, small fissures, dissolved pores and microfractures, there is more contact so temperature variations will be damp and delayed. Amplitude of temperature changes and how quickly it responds to recharge and discharge events can give clues about storage, mixing, and residence time. For example: narrow range of temperature implies well-mixed flow or long-residence time. Even in karst springs with relatively stable mean temperature, careful time-series analysis of temperature can yield insights into aquifer behavior, recharge dynamics and mixing processes [27,28]. Downstream of the karst spring, the water-temperature regime of karst rivers is under the influence of complex heat-exchange process with riverbeds and atmosphere [29]. However, this regime usually is mostly defined by temperature of tributaries and unfortunately anthropogenic inflows, such as wastewater, dam releases and diversions [30,31].

The subject of interest of this paper is the Ruda River located in Croatia, which is characterized by its unique ecological dynamics strongly influenced by anthropogenic factors and karst hydrology. During the last few decades, investigations on this river have been focused on biological and ecological problems. Biological monitoring and physical and chemical spot measurements have been conducted as part of multi-river karst surveys and focused faunal studies. Major studies sampled macroinvertebrates, selected physical and chemical parameters, caddisflies, mayflies and freshwater snails at multiple sites and seasons, producing community and water-quality snapshots rather than long continuous records [32,33,34]. Research on the Ruda River and its catchment highlights ecological sensitivity of karst tributaries and biological responses to flow alteration and local habitat change. Studies show that regulation and direct hydrological alterations can reduce diversity at directly affected sites while lateral springs and karst connectivity can buffer impacts in other reaches [35]. The Ruda is an interesting example of karst river which shows that despite flow regulation the overall species richness and diversity can still be large. The existing research provided good biological baseline data but lacks an integrated systematic assessment of hydrology, which is critical for understanding the karst water system, ensuring safe water supply, protecting ecosystem, and managing sustainably. This paper presents results of hydrological research on the Ruda River based on available data and information. The main aim was the generation of new knowledge about the hydrological functioning of the Ruda River, where the focus was on discharge and water temperature regimes that experienced severe anthropogenic alterations. The methodology used includes classical hydrological and statistical methods as well as methods of time series analysis that, to the best of our knowledge, have not been applied yet in a form like the one presented here.

2. Materials and Methods

2.1. Study Area

The Ruda River is the largest tributary of the Cetina River, which lies in Middle Dalmatia, Croatia (Figure 1). The Ruda Spring is located at the bottom of Kamešnica Mountain at approximate elevation of 298–302.5 m a. s. l. The catchment is set within the Dinaric karst, where conduit and fissure flow in carbonate aquifers produce springs that feed surface streams and create strong groundwater–surface water coupling. It is part of a complex karst hydrological system influenced by groundwater flows, ponors, and subsurface connections [36,37,38,39,40]. It is composed of thick Mesozoic carbonate rocks (limestones and dolomites) with high secondary porosity caused by fracturing and karstification. Dominant lithology is Cretaceous and Jurassic limestones, locally interbedded with dolomitic layers. These rocks are highly permeable due to intense tectonic faulting and dissolution. There are many faults and fractures, forming preferential groundwater flow paths toward the Ruda Spring. Impermeable flysch deposits occur locally, acting as hydrogeological barriers that guide groundwater toward spring zones. Consequently, the aquifer feeding the Ruda Spring is a karst aquifer with conduit and fissure flow. Due to the karstic nature of the terrain, the hydraulic connectivity between surface and subsurface waters is strong and rapid, leading to fast groundwater responses to rainfall events.

Human impacts and hydrological changes are strong, so spring and river function under a substantially altered hydrological regime. The Ruda Spring is vital source of freshwater. Allowed abstraction is 540 l/s or 8,000,000 m³ per year. Actual abstraction is approximately 4,500,000 m³ per year. Before the 1970s, the discharge from spring was controlled by natural karst hydrology, fed by underground aquifers and precipitation cycles. The spring was fed by the water from ponors found in Livanjsko Polje (Bosnia and Herzegovina) and by the infiltrated water from higher karst plateaus of the Kamešnica Mountain. After the Buško Blato reservoir started to work in 1973, significant alterations in natural spring discharge and river flow were seen. The Buško Blato reservoir has an active storage of 786 × 10⁶ m³, with an average level of 716.4 m a. s. l. Water from the Buško Blato reservoir is diverted through the Lipa reservoir and then used for hydropower generation in the Orlovac Hydropower Plant (HP) in Croatia. The Lipa reservoir is situated at 700 m a. s. l., and its capacity is 1.6 million m³. It serves also as a compensation reservoir into which the canal system brings one-day inflows from the Livanjsko Polje. From there, through a 12,100 m long tunnel beneath the Kamešnica Mountain, the water is conveyed to the Orlovac HP. During the rainy period, when the inflow exceeds consumption requirements of the Orlovac HP, the Lipa reservoir fills up, and the surplus water is pumped by a pumping station into the Buško Blato reservoir [41]. Installed capacity of this high-pressure underground hydropower plant is 237 MW, hydraulic head 380 m, and average annual electricity production is approximately 800 GWh. Installed discharge through turbines is 70 m3/s [42]. After being used for power generation, this water is discharged via an open channel into the Ruda River, approximately 1.2 km downstream from the spring. Consequently, Ruda’s flow regime is now mostly artificial, as it receives water primarily from hydropower operations rather than from its natural karst sources. Practically, construction of the Buško Blato and Lipa reservoirs, and operation of the Orlovac HP, have transformed the Ruda River into an engineered discharge outlet, affecting the natural balance and ecological health of the karst river. The artificial regulation decreased natural flow variability (less difference between wet and dry periods), altered temperature and sediment regimes, and reduced ecological connectivity along the river system.

2.2. Available Time Series

Information about precipitations and air temperatures on the catchment of Ruda Spring are obtained from the daily time series collected at the Livno Meteorological Station (MS) located in the Livanjsko Polje and the Rosandići Rain Gauge (RG) situated at the foothills of the Kamešnica Mountain (Figure 1). Information about daily and hourly discharges is obtained from three hydrological stations situated in the wider area of the Ruda Spring: Ruda 1, Ruda Kanal, and Ruda Orlovac (Figure 2). Discharges from the spring are calculated as sum of discharges measured at the hydrological stations Ruda 1 and Ruda Kanal. Cumulative inflows and losses along the segment of river from the Ruda Spring to the Ruda Orlovac Hydrological Station (HS) are represented by differences between measured discharges at the Ruda Orlovac HS and calculated discharges for the Ruda Spring. It should be noted that cumulative inflows and losses are dominantly generated by discharges from the Orlovac HP. However, there are intervals in time series when the discharge at Ruda Orlovac HE is lower than the discharge from Ruda Spring, which denotes the existence of losses at this segment of the river. They are registered during non-operational periods of Orlovac HP. Continuous monitoring of water temperatures is performed only at Ruda Orlovac HS. Basic statistical characteristics of observed air temperatures, precipitations, discharges, and water temperatures are presented in Table 1, including periods of monitoring, analyzed time steps, minimal and maximal registered values, and means.

2.3. Methods

Let $x=\left\{x_1,x_2,\dots ,x_n\right\}$ and $y=\left\{y_1,y_2,\dots ,y_n\right\}$ represent two analyzed time series or two sets of measured data that have means $\mu \left(x\right)$, $\mu \left(y\right)$, standard deviations $\sigma \left(x\right)$, $\sigma \left(y\right)$, and covariance $c\left(x,y\right)$ [43].

Values of Duration Curve (DC) of time series $x$ are estimated by using an empirical cumulative distribution:

$$DC\left(x\right)=\left\{p_i :p_i={\displaystyle \frac{rank\left(x_i\right)}{n+1}}; i=\mathrm{1,2},\dots ,n\right\}$$

(1)

Values of Rescaled Adjusted Partial Sums (RAPS) of time series $x$ are defined as [44,45]:

$$RAPS\left(x\right)= \left\{s_i :s_i={\displaystyle \frac{1}{\sigma \left(x\right)}}\sum _{k=1}^i\left(x_k-\mu \left(x\right)\right); i=\mathrm{1,2},\dots ,n\right\}$$

(2)

Values of Moving Average (MA) filter of time series $x$ for a window size $k$ are calculated by following expression:

$$MA\left(x\right)=\left\{m_i : m_i={\displaystyle \frac{1}{2k}}\sum _{-k}^k{x}_{i+k}; i=k,k+1,\dots ,n-k\right\}$$

(3)

The Pearson correlation coefficient between time series $x$ and $y$ is [43]:

$$r(x,y) = {\displaystyle \frac{c(x,y)}{\sigma \left(x\right)\sigma \left(y\right)}}$$

(4)

Linear Regression (LR) between $x$ and $y$ is defined by:

$$y_i=a{x}_i+b+{\epsilon }_i , i=\mathrm{1,2},\dots ,n,$$

(5)

where $a$ is slope, $b$ intercept, and ${\epsilon }_i$ residuals. Parameters of line, $a$ and $b$, are determined by minimizing the sum of squared residuals. Goodness of fit of LR is estimated by a coefficient of determination $R^2$, where $R=r(x,y)$. Statistical significance of $R$ is evaluated for level $\alpha =0.05$ by using statistics [46]:

$$t={\displaystyle \frac{R\sqrt{n-2}}{\sqrt{1-R^2}}},$$

(6)

that follows a t-distribution with degree of freedom $df=n-2$.

Let $c_k(x,y)$ represent covariance between time series $x=\left\{x_1,x_2,\dots ,x_{n-k}\right\}$, and $y=\left\{y_{1+k},y_{2+k},\dots ,y_n\right\}$, and let $c_k(y,x)$ represent covariance between time series $y=\left\{y_1,y_2,\dots ,y_{n-k}\right\}$ and $x=\left\{x_{1+k},x_{2+k},\dots ,x_n\right\}$, where $k\ge 0$. Values of the Pearson Cross-Correlation Function (CCF) between time series $x$ and $y$ at lag $k$ is defined as (Box and Jenkins, 2008):

$$\begin{gathered} CCF\left(x,y\right)=\left\{r_k : k=0,\pm 1,\pm 2,\dots ,\pm \left|n-k\right|\right\} \\ r_k =\left\{\begin{array}{c}{\displaystyle \frac{c_k(x,y)}{\sigma \left(x\right)\sigma \left(y\right)}} : k\ge 0\\ {\displaystyle \frac{c_{\left|k\right|}(y,x)}{\sigma \left(x\right)\sigma \left(y\right)}} : k<0\end{array}\right. \end{gathered}$$

(7)

3. Results

Analyzed time series are precipitations from Livno MS (P_L), precipitations from Rosandići RG (P_R), air temperatures from Livno MS (T_L), water temperatures from the Ruda Orlovac HS (T_RO), discharges from the Ruda Spring (Q_S), discharges from the Ruda Orlovac HS (Q_RO), and inflows/losses (Q_IL) standing for discharges from the Orlovac HP. Mean, maximal, and minimal values for each year are presented in Figure 3. Homogeneity of time series was controlled by using the Levene test (Brown-Forsythe version) for variances, the independent t-test for means, and the Wilcoxon test for trends [46,47,48]. Results show that annual precipitations P_L and P_R are similar (Figure 3a), but they have a wide range from 761 to 1796 mm. Air temperatures from the Livno MS (Figure 3b) and discharges from the Ruda Spring (Figure 3c) are inhomogeneous. Annual air temperatures T_L have an upward trend, where the increase is approximately 1 °C in 20 years. Annual discharges Q_S have a mean shift in 1988. Namely, the mean annual discharge from Ruda Spring during period 1973-1988 was 14.2 m³/s. During period 1989–2022, it was 12.1 m³/s. The mean shift coincides with the beginning operation of Đale HP situated about 12 km downstream from the spring. Figure 3e presents also results of application of MA filter to annual discharges Q_IL. Results describe inter-annual variations in intensity of operation of the Orlovac HP in the form of five-years means (size of window $k=2$ in Equation 2). Period of the most intensive operation of the Orlovac HP was 1977–1981 with five-years mean of 19.2 m³/s. Period of the lowest intensity was 1984–1988 with five-years mean of only 7.9 m³/s.

RAPS of annual precipitations P_L, and annual discharges Q_RO, Q_S, and Q_IL, are compared in Figure 4. Upward or positive trends in RAPS are periods with values mostly above average. Downward or negative trends are periods with values mostly below average. Sharp change in direction denotes possible regime shift. Maximum or minimum point is potential change-point location. It can be noted that inter-annual variations in annual discharges Q_RO generally have followed variations in annual precipitations P_L during the entire period. On the other hand, two components of Q_RO (Q_S and Q_IL) had different behaviors. They both deviated from annual precipitations P_L, especially during the period before 1988. The change point in year 1988 in RAPS(Q_S) is clear, which confirms inhomogeneity of annual discharges Q_S. Differences in variations between RAPS(Q_S) and other functions are most obvious for period 1983–1988, where only RAPS(Q_S) have positive trends. Concerning annual discharges Q_IL before 1988, RAPS(Q_IL) have positive trends for period 1973–1981, and negative trends after that. Comparing RAPS(Q_IL) in Figure 4 with results of MA filter in Figure 3e, it can be noted that positive trends in RAPS(Q_IL) coincide with the period of the most intensive operation of the Orlovac HP, from 1977 to 1981. After this period, the period of the lowest intensity operation is present, from 1984 to 1988, with negative trends in RAPS(Q_IL). It should be emphasized that negative trends in RAPS(Q_IL) coincide with positive trends in RAPS(Q_S) during this period, which shows that the decrease in annual discharges of the Orlovac HP generated the increase in annual discharges from the Ruda Spring.

DC of daily discharges Q_S for periods 1973–1988 and 1989–2022 are compared in Figure 5a. High-flow segment describes rapid conduit flow, direct recharge response, and flood dynamics. Steeper slope denotes more flashy system. In karst, it is usually very steep, which shows strong hydraulic connection and low attenuation. Middle segment describes mixed flow (conduit + fracture) and transitional storage release. Shape shows degree of internal regulation. Low-flow segment describes matrix drainage, long-term storage and drought behavior. If flat, it shows strong baseflow regulation. If steep, it shows limited matrix storage. The simplest approach to comparing DC is pointwise comparison at characteristic percentiles, P10 for high flows, P50 for mixed flows, and P90 for low flows. It can be noted that corresponding discharges for percentiles P10, P50 and P90 are higher for the first period. In Figure 5b, both DC are presented in log-normal probability plot. Using probability plot enhances visualization of extreme values, makes regime changes more visible, helps find statistical distribution, and simplifies comparison between datasets. Breaks in slope may suggest multiple storage compartments. For example, karst systems often activate conduits above certain recharge thresholds. On probability plot this appears as kink or inflection point or sudden slope increase. This may show overflow conduits, epiphreatic activation, or temporary spring connections. However, DC obtained for period 1989–2022 is rather linear, so almost the entire range of discharges can be approximated by log-normal distribution. Exceptions are only 1% highest values. DC obtained for period 1973–1988 has a more complex that reveals differences in functioning of the karst spring during two periods. The shape of DC for period 1973–1988 at probabilities between 0.1 and 0.99 shows that the karst aquifer of Ruda Spring during this period had more inflows. During both periods baseflow has been stable, which is a characteristic of permanent spring.

RAPS of daily time series including precipitations P_L, discharges Q_S, Q_RO, Q_IL, water temperatures T_RO, and air temperatures T_L are compared in Figure 6. They cover time periods 1996–2021 (P_L, Q_S, Q_RO, Q_IL) or 2009–2023 (T_RO, T_L) depending on comparison requirements and availability of daily data. Since daily time series are autocorrelated, seasonality, trends, and fluctuations are clear in all RAPS. Seasonality dominates in water temperatures and air temperatures, while in precipitations and discharges it is only an additional fluctuation that distorts inter-annual trends. In daily precipitations T_L, three periods can be distinguished, 1996–2009, 2009–2015 and 2015–2021. The first period had below-average precipitations, so RAPS(P_L) has generally negative trends. The second period had the most abundant precipitations, so trends are positive. These trends last until 2015 when they are changed again. The regime of daily discharges Q_S follows the regime of daily precipitations P_L, so inter-annual variations in RAPS(Q_S) and RAPS(P_L) are similar. However, due to retention of rainwater in karst aquifer, sub-annual fluctuations are attenuated, and seasonal variations are amplified, so change points in 2009 and 2015 are evident in RAPS(Q_S). Regime of daily discharges Q_IL is different in such a way that RAPS(Q_IL) is more uniform and without multi-year trends, so a change point does not exist. It is a consequence of Buško Blato reservoir, which enables multi-year regulation of discharges, so daily discharges Q_IL have been generated over the years mostly by electric power demands.

In Figure 6, it can be noted also that variations in RAPS(Q_S), RAPS(Q_IL) and RAPS(Q_RO) have a time delay from variations in RAPS(P_L). The values of these delays are estimated by CCF between two RAPS as the time lag where CCF reaches its maximum. The obtained four CCF are presented in Figure 7. Time delay between RAPS(P_L) and RAPS(Q_S) is an average retention time of rainwater in the karst aquifer of Ruda Spring, and it amounts to 68 days. Time delay between RAPS(P_L) and RAPS(Q_IL) is an average retention time of rainwater that is used by the Orlovac HP. Because of effects of Buško Blato reservoir, this time is larger and amounts 92 days. As expected, the time delay between RAPS(P_L) and RAPS(Q_RO) is between these two values and amounts 82 days. Considering effects of air temperatures on water temperatures, the obtained maximum in CCF(RAPS(T_L), RAPS(T_RO)) has the lowest value of all CCF presented, but it is still statistically significant confirming that air temperatures play a key role. Time delay of 16 days is an average time needed for the heat transfer from the air in Livanjsko Polje to the water in Ruda River by the system of open canals, two reservoirs, and the tunnel of Orlovac HP.

Effects of the Orlovac HP on discharges at the Ruda Orlovac HS are presented in a seasonal daily graph (Figure 8). It shows means as well as minimal and maximal registered discharges for each day of the year (1–365) at the Ruda Spring and the Ruda Orlovac HS during period 1989–2021. These curves show differences in spring and river peaks, summer low flows, autumn response of river and spring to rainfall, and winter baseflow conditions. The regime of river has been primarily defined by operation of HP Orlovac. Maximal daily discharges at the Ruda Orlovac HS were seen during spring months, April and May, up to 108 m³/s. Maximal discharges from Ruda Spring, up to 49.8 m³/s, were seen during autumn months, from October to December, which are most abundant with precipitation. Mean daily discharges have a similar seasonal regime as maximal. Ratio between mean daily discharges MEAN(Q_RO)/MEAN(Q_S) calculated for each day in the year range from 1.34 in September to even 3.00 in Jun. Ratio between maximal discharges MAX(Q_RO)/MAX(Q_S) goes up to 4.88 in July. The regime of minimal discharges has been changed the least, because these values are mostly registered during nonoperation of the Orlovac HP.

The water temperature regime of the river is influenced by the air temperature in Livanjsko Polje and the discharge from Orlovac HP. LR between monthly air temperatures T_L and monthly water temperatures T_RO is presented in Figure 9a. Linear dependence is clear and coefficient of determination $R^2$ is high. However, maximal influence is registered during winter and summer. LR between monthly discharges Q_IL and monthly water temperatures T_RO for February and August are presented in Figure 9b. Values of two intercepts, 9.1 and 15.3 °C, represent average water temperatures for February and August at the Ruda Orlovac HS without effects of discharges from the Orlovac HP.

Concerning seasonal variations in daily air and water temperatures, Figure 10 shows mean values as well as minimal and maximal registered values for each day of the year at the Livno MS and the Ruda Orlovac HS during 2009–2022. Daily air temperatures have clear seasonality. In the coldest month, January, they ranged from -13 to 9 °C with a mean about 1 °C. In the warmest month, July, they were between 12 and 28 °C with a mean of about 22 °C. Irregular fluctuations are evident, especially at minimal values during winter period, from December to February, when air temperatures were between -13 and 1 °C. Two components, seasonal and fluctuations, are observable also in water temperatures, but both are attenuated because of heat capacity of water. During winter months, daily water temperatures ranged from 2 to 10 °C, whereas during summer they were from 12 to 24 °C. Minimal daily water temperatures in the coldest month, January, were between 2 and 4 °C. Maximal daily values were in July and August, between 20 and 24 °C.

Figure 11a and Figure 11b represent basic statistical characteristics of variations in hourly discharges Q_RO and hourly water temperatures T_RO, respectively. Mean values, median, minimal and maximal registered values, as well as percentiles P5 and P95 for each hour of the day (1–24) are presented. The maximal hourly discharge Q_RO during the analyzed period was 122 m³/s. Two peaks in mean discharges, median values, and extreme high discharges P5 are clear, the first at 12 h and the second at 20 h. Since Orlovac HP often works in low operational mode during night hours, discharges at peak hours are more than twice as large as discharges at night hours. For example, ratios between discharges at 20 h and discharges at 3 h are 2.17, 2.46 and 2.12 for mean, median and P5, respectively. Similar two peaks are observable also in water temperatures, at 13 and 21h, most obviously in percentiles P95 and P5 (denoted as T_{RO_95} and T_{RO_5}), showing strong interdependency between variations in discharges and water temperatures. Results reveal also that water temperature in some peak hours during winter can be as low as 0.4 °C, while during summer it can reach up to 25.9 °C.

Similarly, as for Q_RO in Figure 11a, characteristic values P5 and means for each hour in the day have been calculated also for hourly discharges Q_IL. The obtained series are denoted as Q_{IL_5} and Q_{IL_MEAN}. LR between Q_{IL_MEAN} and T_{RO_95}, Q_{IL_MEAN} and T_{RO_5}, Q_{IL_5} and T_{RO_95}, Q_{IL_5} and T_{RO_5}, are presented in Figure 12. Values of coefficient of determination R² are remarkably high, between 0.90 and 0.96. It shows that strong interdependences between hourly discharges from the Orlovac HP and hourly water temperatures at the Ruda Orlovac HS exist from average to extreme values of discharges. Intercepts for discharges Q_{IL_MEAN} (7.2 and 15.9 °C) and Q_{IL_5} (7.5 and 15.2 °C) represent referent values for minimal and maximal water temperatures at the Ruda Orlovac HS without discharges from the Orlovac HP.

4. Discussion

The first step in analyses of time series was identification and explanation of trends and change points. Multiple hydrological and meteorological investigations have demonstrated statistically significant upward trends in air temperatures from meteorological stations in Croatia and Bosnia and Herzegovina [49,50,51,52]. The Livno MS is not an exception, so an upward trend exists. However, there is no statistical evidence that this upward trend in air temperatures in Livanjsko Polje has been reflected in regimes of precipitation on catchment, discharges from the Ruda Spring, or temperatures of the Ruda River.

Effects of hydropower plants on the River Cetina have been elaborated in several studies [37,39,40,42] including effects of the Đale HP on groundwater levels [53]. The mean shift in discharge from Ruda Spring registered in 1988 coincides with the beginning of operation of this hydropower plant. Results show that there was a difference between operation of the Orlovac HP before connection of the Đale HP to hydroelectric power system of the Cetina River and after connection, which most probably affected the functioning of Ruda Spring.

Differences in the functioning of Ruda Spring before and after 1988 were analyzed by using DC of daily discharges. DC is one of the most useful tools in the analysis of karst springs because it condenses long discharge time series into a form that clearly shows variability, regulation capacity, and extremes. DC allows comparison of different karst springs, same spring before and after engineering works, or natural vs. regulated regime [54,55,56]. Probability distributions are sometimes used to parameterize the relationship between discharges and exceeding probabilities, transforming the empirical DC into a continuous analytical function suitable for interpolation, extrapolation, and comparison. When presented in probability plot, DC becomes even more informative for interpreting karst aquifer behavior, so comparison is plainer [57]. In this study, log-normal distribution has been used. The results obtained revealed that DC from the period after 1988 can be approximated by this distribution in almost the entire range of probabilities, from 0.01 to 0.99, which indicates that the Ruda Spring had a relatively consistent functioning without multiple storage compartments (without activations some additional conduits above certain recharge thresholds, epiphreatic activations, temporary spring connections, activations of ponors, etc.). It is not case with the period before 1988, when DC had irregular form that cannot be approximated by any theoretical distribution. Comparison of DC from two periods revealed that DC from the period before 1988 had larger discharges at all segments, low flow, middle flow, and high flow, which implies the existence of continuous inflows during this period. Exceptions are only extreme values, high flows below probability of 0.01 and low flows above 0.99.

Small yet systematic changes over time are often hidden by the variability of data itself but sometimes they enclose particularly essential information [58,59]. RAPS is a statistical method used for detecting these changes in hydrological time series such as precipitation, river discharge, groundwater levels, and water temperature. RAPS visualizes cumulative deviations from the mean value, and thereby highlights change points, data clustering, long-term trends, irregular fluctuations, and periodicities in the time series. Most hydrological papers have applied RAPS to annual precipitation totals or annual mean discharges, because it provides the clearest signal of structural change [44,60,61,62]. Daily values are less proper for classical RAPS trend detection because of strong autocorrelation that produces fluctuations and seasonality. So, RAPS of daily data has been used much less often, for regime change analysis or event-based hydrology [63,64,65]. In this study, RAPS of daily precipitations and discharges were used for comparison of regimes and for estimation of retention times. Results show that inter-annual variations in discharges from Ruda Spring have followed variations in precipitations, while variations in discharges from the Orlovac HP have been more uniform, reflecting effects of water retention in Buško Blato Reservoir. It was noted that time delays between RAPS of precipitations and RAPS of discharges exist, which represent retention times. Average retention times for the entire period were estimated as lags with maximal correlation in CCF between RAPS of precipitations and RAPS of discharges. Because of rainwater retention in the Buško Blato Reservoir, the retention time of rainwater used by the Orlovac HP was 24 days longer than the retention time of rainwater in the karst aquifer of Ruda Spring. The average time needed for heat transfer from the air in Livanjsko Polje to the water in Ruda River was obtained in analogous way by using CCF between RAPS of air temperature and RAPS of water temperature. It shows that changes in the air temperature at the Livno MS have been recorded in the water temperature at the Ruda Orlovac HS after 16 days.

Seasonal effects of the Orlovac HP on the discharge regime of Ruda River were quantified by comparing seasonal daily graphs of the Ruda Spring and the Ruda Orlovac HS. Seasonal daily graphs presented mean values for each day of the year computed over entire observed period, as well as minimal and maximal values registered for each day. They removed inter-annual variability and highlighted the typical seasonal regimes comparing natural vs. regulated discharges. Results show that minimal discharges have increased, total volume of discharges has doubled, but maximal discharges have changed most significantly. They have increased by almost five times in summer months, and they have larger fluctuations due to irregular operation of the Orlovac HP.

The regime of river water temperature is under strong influence of the air temperature in Livanjsko Polje and the discharge from Orlovac HP. LR was used to confirm interdependence, establish linear relationships between variables, control level of statistical significance of relationships, and estimate referent water temperatures (water temperatures without discharges from the Orlovac HP). The established LR between monthly air temperatures at the Livno MS and monthly water temperatures at the Ruda Orlovac HS had coefficient of determination R² about 0.84. LR between monthly discharges from the Orlovac HP and monthly water temperatures at the Ruda Orlovac HS for February and August had values 0.80 and 0.69, respectively. As result of this influence, the water temperature regime of the Ruda River was changed. Depending on discharges, monthly water temperatures may increase during summer, above 20 °C, and decrease during winter bellow 4 °C. On the other hand, two intercepts show that referent water temperatures for February and August were about 9.1 and 15.3 °C, respectively.

As well as discharges, seasonal variations in daily water temperatures at the Ruda Orlovac HS were also analyzed by seasonal daily graphs presenting mean values of temperatures for each day in the year computed over the entire observed period, and minimal and maximal temperatures registered for each day. Results show that minimal and maximal daily water temperatures come from less than 2 °C in January to above 24 °C in July. Mean values are 6 – 8 °C for winter days and 16 – 18 °C for summer days.

In addition to significant seasonal variations, large sub-daily variations in discharges at the Ruda Orlovac HS were also registered. Namely, discharges from the Orlovac HP have two peaks within the day, at 12 and 20 h, where the largest peak discharges are multiple times greater than minimum discharges in night hours. Ratios between discharges in peak hours and night hours go up to almost 30 for median values. Hourly variations in discharges at the Ruda Orlovac HS are typical for a river influenced by peaking hydropower plant. This plant works in a way that rapidly increases or decreases release to meet fluctuations in electricity demand (hydropeaking operational mode). Because of rapid start/stop cycles needed to match energy peak demand, it has highly variable water release with sub-daily changes and large amplitudes. This type of water release alters river systems in ways that can stress ecological and physical processes downstream by altering discharge and water temperature regimes [66,67]. Hydropeaking usually produces thermopeaking – rapid, short-term fluctuations in river water temperature caused by hydropower plant operations [68,69,70,71]. Consequently, strong interdependences between extreme discharges P5 from the Orlovac HP and extreme water temperatures P5 and P95 at the Ruda Orlovac HS were detected, and two peaks within the day in hourly water temperatures at the Ruda Orlovac HS were registered. These two peaks have a time delay of one hour from peaks in discharges. The time delay is a travel time of water from the Buško Blato reservoir to the Orlovac HP. Estimated referent hourly water temperatures were about 7 °C for winter period and 16 °C for summer, which is in accordance with field measurements. Namely, water temperatures measured during 2004 and 2005 about 100 m downstream from the Ruda Spring ranged from 6.8 °C in February to 16.9 °C in August, with average of 11.8 °C [35]. Hourly water temperatures measured at the Ruda Orlovac HS had much wider range, from 0.4 to 25.9 °C, which proves large alterations of temperature regime.

Despite large alterations of discharge and water temperature regime, the Ruda River still is a permanent karst stream with stable physical and chemical conditions. There is no continuous monitoring of the water quality of the Ruda River, but occasional spot measurements showed that physical and chemical characteristics are typical for karst rivers in Dalmatia, strongly influenced by the limestone geology and groundwater input [34]. pH values ranged from 7.45 to 8.55, showing mildly alkaline waters. Electrical conductivity ranged from 282 – 588 µS/cm, showing moderate mineralization. Alkalinity ranged from 140 – 215 mg/L CaCO₃, classifying the water as moderately hard. Dissolved oxygen levels were high (10.4–15.5 mg/L), and oxygen saturation never dropped below 90%.

5. Conclusions

Discharge and water temperature regimes of Ruda River were investigated by using a methodology based on classical hydrological, statistical and time series analysis methods adapted to particularities of study area and available data. RAPS of annual data are used to compare inter-annual variations in precipitations and discharges, and to find change points. Differences in two operational periods of the Orlovac HP, before and after 1988, are clarified by MA filter. Empirical and analytical forms of DC of are used to (1) compare the regime of Ruda Spring before and after 1988, (2) detect the existence of groundwater interchange, and (3) obtain basic information about the duration of diffuse, conduit, and mixed flow at the spring. RAPS of daily data are used to compare regimes of precipitation, discharge, water temperature, and air temperature, where retention times and heat transfer times have been estimated by CCF between two RAPS. LR between extreme discharges P5 from the Orlovac HP and extreme water temperatures P5 and P95 at the Ruda Orlovac HS have been used to (1) confirm the direct interdependence between the discharge from hydropower plant and the river temperature, (2) establish relationships, and (3) estimate referent water temperatures for the Ruda Orlovac HS without discharges from the Orlovac HP. The results obtained advance existing knowledge about the hydrology of Ruda River.

Air temperatures on catchment have statistically significant trend, which is result of global climate change registered at meteorological stations in wider regions. Discharges from Ruda Spring have a mean shift in 1988, which coincides with the beginning of operation of Đale HP. This mean shift is a consequence of differences in operation of Orlovac HP before and after 1988, and it confirms that annual discharges from the spring depend on the operating mode of Orlovac HP.

After 1988, the karst aquifer of Ruda Spring has had persistent storage without significant gains or losses through interchange with neighboring catchments depending on groundwater levels, or gains through allogenic recharge. The Ruda Spring is a permanent karst spring with a stable baseflow and significant seasonal variations depending on recharge conditions. Average retention time of the rainwater is about 68 days.

The discharge regime of Ruda River has been under strong influence of the discharge from Orlovac HP. Discharges from Orlovac HP have been generated over the years mostly by electric power demands. Installed discharge is larger than maximal discharge seen at the spring. The maximal seen hourly discharge at the Ruda Orlovac HS was 122 m³/s, while the maximal discharge from the Ruda Spring during the same period was only 48 m³/s. Sub-annual redistribution of water exists from autumn to spring months. Sub-daily variations with two pronounced peaks at 12 and 20 h are typical for a river influenced by a peaking hydropower plant. Average retention time of the rainwater used by the Orlovac HP is about 92 days.

The water temperature regime of the Ruda River has been substantially altered. Depending on discharge, monthly water temperature ranges from bellow 4 °C during winter to above 20 °C during summer. In addition to seasonal variations, hydropeaking produced thermopeaking so there are two peaks in sub-daily variations of water temperature. Since the water from Buško Blato reservoir needs time to reach the Orlovac HP, these two peaks have a time delay of about one hour from peaks in discharge. Hourly temperatures seen at the Ruda Orlovac HS ranged from 0.4 to 25.9 °C. Without discharges from the Orlovac HP, water temperatures would range from about 7 °C to about 16 °C. Changes in the air temperature in Livanjsko Polje are observable in the river water temperature after 16 days.