Evaluation and Quantiﬁcation of the E ﬀ ects of Climate and Vegetation Cover Change on Karst Water Sources: Case Studies of Two Springs in South-Western Slovenia

: Karst aquifers hold important water resources such as regional water budgets and freshwater supply. Due to highly dynamic hydrological processes in comparison to other less permeable systems, they are particularly susceptible to environmental changes. However, little research directly characterizes the impacts of climate and vegetation cover changes on karst water sources. In this study, we aimed to evaluate individual long-term impacts and impacts of sudden large-scale forest disturbances on changes in groundwater recharge and in spring discharge. The work is based on temporal analysis of forest cover and a comparison of trend analysis of hydro-meteorological parameters. To investigate and evaluate vegetation cover change impacts on groundwater recharge, we used a soil water balance model and compared estimated actual daily values of e ﬀ ective precipitation to its ﬁctional estimation disregarding the vegetation cover change. The applied methodology enabled quantiﬁcation of the impacts of climate and vegetation cover change on selected karst water sources. The study suggests that the vegetation cover can have a signiﬁcant impact on the spring recharge. Large-scale disturbances that occurred in a short-term mitigated the e ﬀ ects expected from the trend analysis of hydro-meteorological parameters. In the long-term, in addition to climate changes, the multi-decadal natural vegetation overgrowth signiﬁcantly contributed to the reduction in the spring’s discharge values, especially in the warm season when water demand is higher. Therefore, the results are of key importance for developing proper water management and environmental policies.


Introduction
The importance of aquifers with karst porosity, which already supply around a quarter of the global needs for drinking water, is gradually increasing [1]. This proportion even exceeds half of the demand of individual countries in the areas of the Dinaric Karst [2], the largest contiguous karst region in Europe [3], and in the Alps. These areas are mostly covered by forest and are important for its rich and unique (underground) ecosystems [4,5].
Due to the dissolution of the host rock, the karst surface is very permeable, the soil cover is usually clayey, very thin or even absent. Characteristic features are the immediate infiltration data of high temporal resolution; (3) aimed to assess the direct impact of changes in climate and vegetation cover on karst water sources. More specifically, the aim of the study was to compare longand short-term changes in groundwater recharge (i.e., effective precipitation) and in spring discharge, thus assessing the response of a karst aquifer to climate and vegetation cover changes.

Study Areas
Two binary karst catchments in SW Slovenia have been selected as representative cases for the study (Figure 1, Table 1). The catchment area of the Unica karst springs enabled the study of the effects of large-scale forest disturbances (i.e., ice breakage, bark beetle infestation, windthrow) in the period 2014-2018 on the discharge conditions of the springs. The catchment area of the Rižana karst spring has been influenced by climate and land use changes over the last six decades and enabled studies of their long-term effects. The catchment of the Unica springs extends over 800 km 2 and embraces three distinctive subcatchments [55][56][57]; karst poljes of the Notranjska, the Pivka River basin at altitudes between 500 and 700 m, and the Javornik-Snežnik karst plateau, which reaches an altitude of 1800 m. In the first two subcatchments, dominant rock types are dolomite or flysch covered with eutric brown soils, which induce surface runoff that is part of the catchment area of the investigated spring. On the karst plateau, highly karstified limestone with very thin rendzinas or brown calcareous soil prevails, which conditions the direct infiltration of precipitation underground [58]. The catchment belongs to the moderate continental climate and is mostly covered with Abieti-Fagetum dinaricum forests. A greater concentration of settlements, industry and agricultural activities is at the bottom of the karst poljes and in the Pivka basin. The flow rates of the Unica River in the period 1989-2018 were between 1 and 90 m 3 /s, while the mean flow rate was 21 m 3 /s [59]. The average annual precipitation in the same period was estimated at 1505 mm, the average annual potential evapotranspiration (ETP) at 738 mm, the average annual runoff at 767 mm and the average annual runoff ratio at 0.51 [60].
The catchment of the Rižana spring is estimated at 247 km 2 [61,62]. The predominant part cover highly karstified limestones and rendzina soil. The area is mainly forested with Orno-Quercetum petraeae-pubescentis. The spring is also recharged by the sinking rivers at the southern margin of the flysch Brkini Hills, where distric brown soils prevail. The location of the spring is at an altitude of 70 m while the catchment lies at between 500 and 1000 m altitudes. The area mainly belongs to the moderate Mediterranean climate. The flow rates of the Rižana River in the period 1989-2018 were between 70 L/s and 63 m 3 /s, while the mean flow rate was 3.5 m 3 /s [59]. The average annual amount of precipitation in the same period was 1484 mm, the average annual ETP was 892 mm, the average annual runoff was 592 mm and the average annual runoff ratio 0.4 [60].

Input Data
Meteorological and hydrological data were obtained from the archive of the Slovenian Environment Agency's website [59,60]. For analysis of the Unica springs catchment, data on total monthly and annual precipitation and ETP from the Postojna meteorological station and data on monthly and annual characteristic discharges from the gauging station Hasberg on the Unica River for the period 1962-2018 (57 years) were used. For analysis of the Rižana spring catchment, data on total monthly and annual precipitation from the Podgrad precipitation station and data on monthly and annual characteristic discharges from the gauging station Kubed II on the Rižana River for the period 1965-2018 (54 years) were used. Data on total monthly and annual ETP for the Podgrad meteorological station was calculated as a mean value regarding data from the Postojna and Portorož-Airport meteorological stations. Podgrad precipitation station lies outside the Rižana spring catchment, yet it is the most representative one for the catchment. ETP is not a directly measured meteorological variable, but is calculated by the widely used Penman-Monteith equation recommended by FAO, including daily data on maximum and minimum air temperature, relative humidity, global net radiation and wind speed [63]. For the assessment of effective precipitation (P ef ) in the Rižana and Unica catchments in the period 2014-2018, daily data on precipitation, snow cover depth, air temperature, average wind speed, relative humidity, and sunshine hours from the Postojna, Portorož-Airport and Kozina meteorological stations were used. For rain events within individual days, an average rain intensity was defined on the basis of precipitation data at 30-min intervals from the stations Postojna and Portorož-Airport.
Vegetation cover data were obtained from the historical digital orthophotographs from 1957 and recent national database on land use from 2013 and from 2018 [64]. Data on growing stock for the period 2014-2018 were provided by Slovenia Forest Service (unpublished data) and data on soil characteristics were extracted from soil maps [65].

Statistical Analysis
In the study, significant trends of change in meteorological and hydrological variables were determined using the non-parametric Mann-Kendall method [66,67], as it is distribution-free, more robust to outliers, and has a higher power than many other commonly used tests [68]. For the true slope of an existing trend (as change per year) the non-parametric Sen's slope method (the Theil-Sen estimator) was used, which is one of the most commonly used tests for estimating linear time trends [69][70][71][72][73]. Slope estimator (β) as a measure of trend magnitude is calculated as the median value over all pairs of points in the time series [68]. Since this approach is more robust to outliers, compared to the standard linear regression, it has been widely used in identifying the slope of the trend lines of meteorological and hydrological time series (e.g., [12,[74][75][76][77][78][79][80][81][82][83][84]). In the present study, we consider results of the trend analysis as statistically significant in cases where p-value exceeds 0.1. One exception refers to the mean annual discharge of the Unica River in the period 1962-2018, where we consider the p-value of 0.1084 as statistically significant, too. Further explanation can be found in the Section 3.1 and a more detailed explanation in the Supplementary Materials.

Vegetation Cover Changes Analysis
In the Unica catchment, the primary focus were effects of large-scale forest disturbances such as ice breakage, bark beetle infestation and windthrow that occurred in the period 2014-2018. For these purposes, a national database on land use for the years 2013 and 2018 was used. However, according to these maps, the forest area in the catchment area of Unica springs remained practically the same in the period 2014-2018. Therefore, additional information was gathered from the annual data on the growing stock provided by the Slovenia Forest Service. Each forest section with growing stock change was considered as damaged or lost forest area on an annual basis. The gained information was used to evaluate annual reduction in forest density in the observed period. This decrease was taken into account in the assessment of the forest cover in the Unica springs catchment area.
In the Rižana catchment, interest was shown in the long-term effects of changes in vegetation cover to karst spring discharge conditions. Therefore, different categories of land use cover were obtained from digital orthophotographs from 1957 and compared to the recent national database on land use from 2018 [64].
Changes in vegetation cover were determined with the help of ArcGIS 10.7.1 Software (Esri, GDi GISDATA, Ljubljana, Slovenia). The obtained information for both study areas enabled us to distinguish between areas with and without forest.

Assessment of the Influence of Changes in Vegetation Cover on the Effective Precipitation
As the main objective of the study was to make a general assessment of whether major changes in vegetation cover cause perceptible changes in the recharge of karst aquifers, and consequently in the discharge of drinking water sources, several existing tools for this assessment have been considered. However, a semi-distributed model VarKarst [85], which simulates karst potential recharge, cannot separate the effects of climate and vegetation cover changes. According to the authors of the LuKARS model [86,87], this semi-distributed model is more practical for catchments of smaller size, and therefore not very suitable for application in catchments of several 100 km 2 . For the purposes of this study, the V2Karst model [40], which is specified for integrated vegetation-recharge modelling in karst, would be suitable.
Due to the large extent of the study areas, their complex hydrogeological settings and the limited data set available, a soil water balance model was used. This decision also resulted from the fact that the model was first developed and its performance tested in relation to the data obtained in the study of the catchment area of the Vipava karst spring in western Slovenia, which has very similar meteorological, hydrogeological and vegetation characteristics to the two study areas [88,89]. The P ef , which represents the amount of precipitation that actually penetrates from the soil layer into the rock layer, was assessed as an average value for the entire catchment area at daily intervals. Only two basic vegetation types were considered: forest and area without forest. This is acceptable as more than 2/3 of the study areas are covered by forest and the research was focused on the effects of forest cover change. The main equation used was The source of water is precipitation, reduced by the amount of precipitation, collected by the vegetation cover (P g ). Under favorable climatic conditions, the additional quantities of water are provided by the snowmelt (M). The actual evapotranspiration (ETR), through which the water is returned to the atmosphere, thus also denotes the amount of water consumption. The relationship between absorption and consumption is influenced by the water storage capacity of the soil, which in turn determines the amount of water that can be stored in it. The change in water storage in the soil is defined by the parameter ∆S.
The precipitation data were corrected for errors due to wind and wetting losses [90]. Interception is defined as the process by which precipitation falls on vegetative surfaces, while loss due to this interception is the proportion of water that returns to the atmosphere by evaporation [91]. For the assessment, we used the conceptual Rutter model [92], which calculates a running water balance of the canopy and tree trunks and separates the areas with and without forest. Experimental data on parameters describing the properties of the canopy and the tree trunks were not available. Therefore, they were introduced into the model as calibration parameters with different values for the foliated and leafless periods. The empirical equations introduced by the US Army Corps of Engineers [93] were used to calculate the amount of snowmelt. They are presented in different forms depending on the precipitation and vegetation conditions.
The modified Penman equation [94] was used to assess ETP. In this way, it was possible to estimate the amount of precipitation reaching the ground, the amount of melted snow and ETP.
On the other hand, the evaluation of the ETR and the changes in water storage in the soil depend on the hydrological properties of the soil. The basic soil properties were obtained from the soil maps of the study areas [65], but the soil hydrological parameters were subsequently corrected in the calibration process. The calibration process was based on the general equation of the water budget: water input = water output ± change in storage, where the water input is the volume of water produced by P ef in the recharge area and the water output is the volume of water discharged by the spring within a selected time interval. A certain error has been made in assuming that the change in storage is negligible. However, this simplification was considered acceptable, since the observed interval of 5 years starts and ends with similar hydrological conditions at the end of a recession period with very similar slopes of the recession curve and values of the minimum discharge ( Figure S1 in Supplementary Materials). The delineation of the two catchment areas is based on many hydrogeological studies, including numerous tracer tests. Since the directions and characteristics of groundwater flow are well known, the possibility of water balance errors due to overlooked inter-basin groundwater flows in both large catchments is significantly reduced. However, such errors cannot be completely avoided due to the complex structure and functioning of karst aquifers.
The model, which combines all the above-described calculations, was set in a Microsoft Excel file [88]. A more detailed description of the parameters and equations used can be found in the Supplementary Materials. First, the model was used in the catchment of the Unica springs for the assessment of the impact of large-scale disturbances on effective precipitation. For the period 2014-2018, it was calibrated (more details in Supplementary Materials) based on the actual meteorological, hydrological and land cover data (gradual decrease in forest cover due to large-scale disturbances). The model was then used with the assumption that the large-scale disturbances did not occur. In this way, the calculated notional amount of total P ef was compared with the actual amount.
In the catchment area of the Rižana spring, the model was used to assess the effects of natural vegetation overgrowth on P ef . For the period 2014-2018, it was calibrated (more details in Supplementary Materials) on the basis of the actual meteorological, hydrological and land cover data. Then, it was used under the assumption that natural overgrowth did not occur and the share of forest remained the same as in 1957. In this way, the calculated notional amount of total P ef was compared with the actual amount.

Trend Analysis of the Hydro-Meteorological Parameters
In case of the Unica River, long-term trend analysis in the period 1962-2018 confirmed statistically significant decreasing trend in annual P (28 mm per decade; p < 0.05), increasing trend in annual ETP (28 mm per decade; p < 0.001) in its catchment.
In the 57-year measuring period (1962-2018; Qmean = 21.967 m 3 /s) of the mean annual discharge of the Unica River at the Hasberg gauging station (Slovenia), in 2014, a record mean yearly discharge of 40.76 m 3 /s was measured. This is the only value in the entire observation period that is greater than the sum of (1) the mean value in the period (21.967 m 3 /s) and (2)  2014 in the trend calculations (the methodology we used allows this), the p-value of the mean annual discharge of the Unica River was lower than 0.05.
Overall, we decided to consider the mean annual discharge trend of the Unica River as statistically significant. It shows a decreasing trend of 820 L/s per decade (Table 2, Figure 2). Statistically significant are also negative trends of the absolute minimum annual discharge (170 L/s per decade; p < 0.05).
In case of the Rižana River, long-term trend analysis in the period 1965-2018 confirmed a statistically significant increasing trend in annual ETP of 35 mm per decade (p < 0.001) in its catchment, while the decreasing trend in annual P (39 mm per decade) is not statistically significant (Table 2; Figure 3). Long-term changes in values of P and ETP are reflected in the reduced discharge of the Rižana River, for which a statistically significant decreasing trend in mean annual discharge of 380 L/s per decade (p < 0.001) is characteristic (Figure 3). The decreasing trends of the absolute minimum annual discharge (30 L/s per decade; p < 0.001) and a maximum annual discharge (2.75 m 3 /s per decade; p < 0.05) are statistically significant.

Evaluation of Changes in Vegetation Cover and Their Influence on Effective Precipitation
Although according to the national database on land use the area of forest in the catchment of the Unica springs in period 2014-2018 remained practically the same, the forestry data show that growing stock decreased for 16.77% due to the large-scale forest disturbances. On the basis of this information, it was calculated that in the assessment of P ef , the share of forests should be gradually reduced from 77.1% before 2014 to 64.2% by the end of the 2014-2018 period. These values were implemented in the soil water balance model in Microsoft Excel and the model was calibrated (calibration parameters: soil water storage capacity and canopy properties) to correspond to the mean discharge of the Unica River in the same period (21.24 m 3 /s). In this way, the total P ef for the period 2014-2018 was estimated at 4875 mm (annual average 975 mm). Assuming that the large-scale disturbances did not occur, in the model the share of forest was set to 77.1% over the entire period. This would reduce the total P ef in the years 2014-2018 to 4643 mm (annual average 929 mm).
A comparison of the monthly values of P ef calculated with the two models ( Figure 4) shows that the ratio ranges between 1 and 1.6 and depends on the amount of precipitation and the season. As expected, in dry periods with active vegetation, the increase of P ef due to the decrease in the share of forests is evidently higher.
In the catchment of the Rižana spring, at present more than 2/3 of the total area is covered with forests. Only about 60 years ago, this share was slightly above 1/3. To evaluate the influence of such intensive natural vegetation overgrowth on the recharge of the karst aquifer, the same calculation mode was used as for the Unica karst spring. The actual share of forest (67.6% of the total catchment area) was taken into account and the model was calibrated to correspond to the mean discharge of the Rižana River in the period 2014-2018 (4.12 m 3 /s). In this way, the total P ef for the period 2014-2018 was estimated to be 3579 mm (annual average 716 mm). Assuming that natural overgrowth did not occur and the share of forest remains the same as in 1957 (37.4%), this amount would increase by 20% to 4293 mm (annual average 857 mm).
Comparison of the monthly values of P ef calculated with the two models ( Figure 5) shows that the ratio ranges between 0.5 and 1. The ratio depends on the amount of precipitation and the season. As expected, in dry periods with active vegetation the natural vegetation overgrowth has higher influence and significantly reduces the amount of P ef . Table 2. Values of Theil-Sen slope estimator (β) for selected meteorological and hydrological variables with labelled statistical significance. Units of precipitation (P) and ETP are in mm and of mean (Qmean), minimum (Qmin) and maximum (Qmax) discharge are in m 3 /s. 2014 in the trend calculations (the methodology we used allows this), the p-value of the mean annual discharge of the Unica River was lower than 0.05. Overall, we decided to consider the mean annual discharge trend of the Unica River as statistically significant. It shows a decreasing trend of 820 L/s per decade (Table 2, Figure 2). Statistically significant are also negative trends of the absolute minimum annual discharge (170 L/s per decade; p < 0.05). In case of the Rižana River, long-term trend analysis in the period 1965-2018 confirmed a statistically significant increasing trend in annual ETP of 35 mm per decade (p < 0.001) in its catchment, while the decreasing trend in annual P (39 mm per decade) is not statistically significant (Table 2; Figure 3). Long-term changes in values of P and ETP are reflected in the reduced discharge of the Rižana River, for which a statistically significant decreasing trend in mean annual discharge of 380 L/s per decade (p < 0.001) is characteristic (Figure 3). The decreasing trends of the absolute minimum annual discharge (30 L/s per decade; p < 0.001) and a maximum annual discharge (2.75 m 3 /s per decade; p < 0.05) are statistically significant.

Evaluation of Changes in Vegetation Cover and Their Influence on Effective Precipitation
Although according to the national database on land use the area of forest in the catchment of the Unica springs in period 2014-2018 remained practically the same, the forestry data show that growing stock decreased for 16.77% due to the large-scale forest disturbances. On the basis of this information, it was calculated that in the assessment of Pef, the share of forests should be gradually reduced from 77.1% before 2014 to 64.2% by the end of the 2014-2018 period. These values were implemented in the soil water balance model in Microsoft Excel and the model was calibrated (calibration parameters: soil water storage capacity and canopy properties) to correspond to the mean discharge of the Unica River in the same period (21.24 m 3 /s). In this way, the total Pef for the period

Discussion
Both the Unica and the Rižana springs show a negative long-term trend in mean annual discharge in the studied periods. However, results of the long-term trend analysis show that the magnitude of a decrease in runoff differs between the studied springs. Namely, the relative decrease in mean annual discharge is greater for the Rižana spring. Since an apparent difference in long-term vegetation cover change in the catchments of the studied spring exists, both climate and vegetation factors have been considered for the interpretation of a long-term decrease of mean annual discharges of the springs. The effects of large-scale forest disturbances have been studied and discussed in the catchment of the Unica springs. Mean monthly discharges of the Unica River in the cold span of the year show an increasing trend that is most probably the result of an increase in air temperatures in winter months in the region (Table 2; Figure 2). Higher temperatures mean less snow precipitation and a shorter snow residence time, which is reflected in increased runoff during the cold period of the year, when, in addition, vegetation is not active; something that has already been observed in catchments with similar relief and climate characteristics elsewhere [83,[95][96][97][98]. The effect of the described phenomenon on runoff in the Unica catchment is not significantly hindered by the increased ETP. Namely, monthly ETP trends from December to March show a statistically significant increase.

The Unica Case Study
Long-term changes in the vegetation cover in the catchment of the Unica springs are relatively small for the entire study period 1962-2018. As these only account for a few percentage points, changes in the vegetation cover were not attributed any major impact on the springs. The impact of climatological changes (in particular, the increase in ETP) was instead considered more influential in the overall decrease in the Unica River discharge, although some studies indicate that there are significant seasonal variations in ETP, especially in temperate climatic zones where deciduous trees predominate and are not often exposed to drought conditions [40,99,100].
More significant are the vegetation cover changes caused by ice breakage, bark beetle infestation and windthrow, which occurred in the years 2014-2018, with a gradual decrease in the proportion of forest from 77.1% to 64.2%. We therefore hypothesized consequent changes in the rate of precipitation interception and transpiration. A higher infiltration rates of precipitation into unsaturated zone and subsequent higher recharge of groundwater was expected.
To assess influence of altered forest cover on P ef , two different scenarios were tested in a model of the soil water balance. Results show that the large-scale disturbances caused an increase in P ef of 5%; taking into account the catchment size, this would mean an increase in mean discharge of 1 m 3 /s. We explain this with the fact that reduced forest vegetation caused diminution of canopy interception, water consumption and soil evaporation (Figure 6), which are dominant components of evapotranspiration [99,101].
Changes in evapotranspiration processes and consequent groundwater recharge due to altered forest cover have previously been shown by other studies as well [100,102]. Although our study provided a rough estimate that cannot take into account the exact natural conditions, the results are consistent with the comparison discussed hereafter. For the 5-year period 2014-2018, compared to the 52-year period 1962-2013, an increase in annual precipitation of 2.7% and annual evapotranspiration ETP of 14.8% and a decrease in annual runoff, calculated as the direct difference between precipitation and ETP, of 6.9% was calculated. Changes in the values of the above-mentioned water balance elements should be reflected in a significant decrease in the mean discharge of the Unica River. However, in the period 2014-2018 its mean discharge decreased by only 0.4% compared to the period 1962-2013. The results obtained support the assumption that the decrease in forest density due to large-scale disturbances causes an increase in P ef and a relative increase in the discharge of the Unica River. Furthermore, our findings are consistent with virtual experiment conducted by Sarrazin et al. [40], which is focused on the impact of land cover change.
Water 2020, 12, x FOR PEER REVIEW 13 of 20 forest density due to large-scale disturbances causes an increase in Pef and a relative increase in the discharge of the Unica River. Furthermore, our findings are consistent with virtual experiment conducted by Sarrazin et al. [40], which is focused on the impact of land cover change.

The Rižana Case Study
From the beginning to the end of the observed period 1965-2018, the mean annual discharge of the Rižana River dropped from 4.75 to 2.74 m 3 /s, or by 51% according to the trend function. For the spring catchment, using the trend function, we determined a 13.7% drop in annual precipitation from 1625.96 to 1416.08 mm in the same period. Since the negative precipitation trend is not statistically significant, we cannot consider it as an important factor of the drop in mean annual discharge of the Rižana River in the observed period, as was already indicated for the Unica River. Furthermore, unlike the Unica River, where we confirmed an increasing trend in mean discharges in winter months, no such phenomenon was determined for the Rižana River. In comparison to the Unica springs catchment, the Rižana spring catchment lies on lower altitudes with less snow precipitation, hence the effect of decreasing snow precipitation and a consequently shorter snow residence time due to a general increase in air temperatures in the observed period was not so evident [25]. Concerning water availability for the drinking water supply of the region in the warm season, when the water demand is the highest, statistically significant negative monthly trends of Qmean and Qmin of the Rižana spring from April to September are discouraging (Table 2; Figure 3).
Although Ravbar et al. [25] gave the effects of climatological changes as the reason for a decrease in mean annual discharge of the spring, this study seeks additional causes. The observed significant drop in the mean annual discharge cannot be simply explained by the effects of climate change with increased annual ETP in its catchment that, in the observed period, increased from 756.65 to 942.80 mm or for 22%, which is exactly the same ratio as we calculated for the Unica catchment.
A possible explanation for such a huge decrease in mean annual discharge of the Rižana River is a distinctive change in vegetation cover in its catchment in the last six decades. The percentage of forests increased from 37.4% in 1957 to 67.6% in 2018, which could have a significant impact on changing the recharge characteristics of the entire catchment. As shown in a study by Huxman et al. [103], woody plant encroachment leads to increased plant transpiration and simultaneously to reduced availability of soil water. Similar anticipations have been reported by previous studies undertaken in karst areas [99,104].
Analysis of evaluation of Pef in the period 2014-2018 for the Rižana spring, assuming that natural overgrowth did not occur and the share of forest remained the same as in 1957, shows that the amount of Pef would be 20% higher and supports the assumption that forest cover expansion leads to a reduction in groundwater recharge. The findings also show that without significant natural

The Rižana Case Study
From the beginning to the end of the observed period 1965-2018, the mean annual discharge of the Rižana River dropped from 4.75 to 2.74 m 3 /s, or by 51% according to the trend function. For the spring catchment, using the trend function, we determined a 13.7% drop in annual precipitation from 1625.96 to 1416.08 mm in the same period. Since the negative precipitation trend is not statistically significant, we cannot consider it as an important factor of the drop in mean annual discharge of the Rižana River in the observed period, as was already indicated for the Unica River. Furthermore, unlike the Unica River, where we confirmed an increasing trend in mean discharges in winter months, no such phenomenon was determined for the Rižana River. In comparison to the Unica springs catchment, the Rižana spring catchment lies on lower altitudes with less snow precipitation, hence the effect of decreasing snow precipitation and a consequently shorter snow residence time due to a general increase in air temperatures in the observed period was not so evident [25]. Concerning water availability for the drinking water supply of the region in the warm season, when the water demand is the highest, statistically significant negative monthly trends of Qmean and Qmin of the Rižana spring from April to September are discouraging (Table 2; Figure 3).
Although Ravbar et al. [25] gave the effects of climatological changes as the reason for a decrease in mean annual discharge of the spring, this study seeks additional causes. The observed significant drop in the mean annual discharge cannot be simply explained by the effects of climate change with increased annual ETP in its catchment that, in the observed period, increased from 756.65 to 942.80 mm or for 22%, which is exactly the same ratio as we calculated for the Unica catchment.
A possible explanation for such a huge decrease in mean annual discharge of the Rižana River is a distinctive change in vegetation cover in its catchment in the last six decades. The percentage of forests increased from 37.4% in 1957 to 67.6% in 2018, which could have a significant impact on changing the recharge characteristics of the entire catchment. As shown in a study by Huxman et al. [103], woody plant encroachment leads to increased plant transpiration and simultaneously to reduced availability of soil water. Similar anticipations have been reported by previous studies undertaken in karst areas [99,104].
Analysis of evaluation of P ef in the period 2014-2018 for the Rižana spring, assuming that natural overgrowth did not occur and the share of forest remained the same as in 1957, shows that the amount of P ef would be 20% higher and supports the assumption that forest cover expansion leads to a reduction in groundwater recharge. The findings also show that without significant natural vegetation overgrowth the average discharge of the spring would increase from 4.12 to 4.95 m 3 /s. This value is in accordance with the calculated mean annual discharge according to the trend function of the spring for the year 1965 (4.75 m 3 /s). If the comparison is focused only on the summer months (June, July, August), the difference in P ef is much higher, reaching up to 49%.
The statistically significant decrease in mean annual discharge of the Rižana River in the period 1965-2018 can thus be considered as a result of interacting impacts of long-term climate and vegetation cover changes (Figure 7), which is in line with virtual simulations made by Sarrazin et al. [40]. This value is in accordance with the calculated mean annual discharge according to the trend function of the spring for the year 1965 (4.75 m 3 /s). If the comparison is focused only on the summer months (June, July, August), the difference in Pef is much higher, reaching up to 49%. The statistically significant decrease in mean annual discharge of the Rižana River in the period 1965-2018 can thus be considered as a result of interacting impacts of long-term climate and vegetation cover changes (Figure 7), which is in line with virtual simulations made by Sarrazin et al. [40].

Conclusions
Our research is a comprehensive study that evaluates the changes within the entire catchment area in order to investigate individual influences, i.e., of climate and vegetation cover changes, on groundwater recharge in karst regions. It suggests that the vegetation cover change in a catchment of a karst spring can, in addition to climate changes, have a significant impact on the spring hydrology over a short and long timespan. Based on observational data, a further finding of this study is that the methodology applied allowed the quantification of impacts of these changes on karst water sources.
The case of the Unica springs catchment enabled the study of the effects of changes in vegetation cover to karst spring discharge conditions due to large-scale forest disturbances in the period 2014-

Conclusions
Our research is a comprehensive study that evaluates the changes within the entire catchment area in order to investigate individual influences, i.e., of climate and vegetation cover changes, on groundwater recharge in karst regions. It suggests that the vegetation cover change in a catchment of a karst spring can, in addition to climate changes, have a significant impact on the spring hydrology over a short and long timespan. Based on observational data, a further finding of this study is that the methodology applied allowed the quantification of impacts of these changes on karst water sources.
The case of the Unica springs catchment enabled the study of the effects of changes in vegetation cover to karst spring discharge conditions due to large-scale forest disturbances in the period 2014-2018. The findings show a resulting increase in P ef by 5% (due to less canopy interception and evapotranspiration) and a relative increase in the discharge of the Unica River, which were due to changes in the water balance elements, expected to decrease more significantly. The effects of large-scale forest disturbances thus mitigated the effects expected from the trend analysis of hydro-meteorological parameters.
The catchment of the Rižana spring has been impacted by climate and land use changes in the past six decades and enabled studies of their long-term impacts. Comparing the effects of climate change, reflected in the long-term increase in ETP, and of change in percentage of forest, resulting in a long-term decrease in P ef , the study points to interacting impacts of climate and vegetation cover changes. Without significant natural vegetation overgrowth, the mean discharge values of Rižana would be 4.95 m 3 /s, which is 0.83 m 3 /s higher than the actual average value. Even more significantly, up to 49% higher P ef and consequently higher discharges would be expected in summer months (June to August), when water scarcity is a serious problem under current conditions.
Although the study is site-specific, the selected karst catchments were nevertheless recognized for the following contributions: (1) well studied complex binary karst aquifers, (2) containing a long history of both climatological and hydrological datasets of high quality and (3) exhibiting strong and clear long-term vegetation cover changes or large-scale forest disturbances. Therefore, the results of this research are transferable to other similar studies aiming to assess the direct impact of various environmental changes on karst water sources. The study is also useful for the further evaluation of ecohydrological processes and the improvement in environmental policies, as it was able to distinguish between effects of climate and land cover changes on hydrological processes in karst. Furthermore, the study has a direct relevance to numerous land management challenges, due to issues of water availability to water-dependent economic sectors in particular.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-4441/12/11/3087/s1, Figure S1: Cumulative P ef for the Unica study area in the period 2014-2018 calculated in the calibration process using various combinations of calibration parameters, Figure S2: Cumulative P ef for the Rižana study area in the period 2014-2018 calculated in the calibration process using various combinations of calibration parameters, Table S1: Measured parameters, Table S2: Parameters taken from published tables and previous applications of the model in areas with similar hydrogeological, vegetation and climate conditions, Table S3: Calibration parameters, Table S4: Calibration process.
Author Contributions: All authors participated significantly and equally in the conceptual and methodological preparation of the paper, in data collection, their analysis and interpretation, as well as contributed to graphical production. Trend analyses were made by G.K., modelling was conducted by M.P. Project administration and funding acquisition was done by N.R. All authors have read and agreed to the published version of the manuscript.