# A New Index to Assess the Effect of Climate Change on Karst Spring Flow Rate

Ahmad Behrouj Peely, Zargham Mohammadi, Vianney Sivelle, David Labat, Mostafa Naderi

► To cite this version:

Ahmad Behrouj Peely, Zargham Mohammadi, Vianney Sivelle, David Labat, Mostafa Naderi. A New Index to Assess the Effect of Climate Change on Karst Spring Flow Rate. Sustainability, 2024, 16 (3), pp.1326. 10.3390/su16031326. hal-04818777

HAL Id: hal-04818777

https://hal.science/hal-04818777v1

Submitted on 5 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# A New Index to Assess the Effect of Climate Change on Karst Spring Flow Rate

Ahmad Behrouj Peely <fnref n="1" />,* , Zargham Mohammadi <fnref n="1" /> , Vianney Sivelle <fnref n="2" />, David Labat <fnref n="3" /> and Mostafa Naderi <fnref n="4" />

**Abstract:** Karstic aquifers, because of their conduit system, are susceptible to climate change. Ten karst springs in the Zagros region were selected to investigate the impact of climate change under three CMIP6 scenarios: SSP1-1.9, SSP2-4.5, and SSP5-8.5. This study was conducted in three steps: downscaling climate projection, analyzing spring discharge time series, and introducing a new index to assess the impact of climate change on spring flow rate. Applying LARS-WG6, precipitation was downscaled at 14 stations in the study area. Moreover, time series and trend analysis showed that the selected springs have experienced a decrease in their flow rate. Assuming the covariance function between precipitation and spring discharge is constant, new indices (i.e., *I*<sub>Qd</sub>, *I*<sub>dQd</sub>, and *I*<sub>cc</sub>) were introduced to highlight the effect of climate change according to the three scenarios. *dQd* is the variability of spring discharge from past to future, *I*<sub>dQd</sub> is spring discharge variability over the historical data, and *I*<sub>cc</sub> is the effect of precipitation and spring discharge change together. *I*<sub>cc</sub> has a range from -0.25 to 0.25 below and above, which is indicative that two extreme conditions including the spring dryness and overflow are in effect, respectively. The main results revealed that the degree of impact at each spring is a function of climate change scenarios and hydrogeological characteristics of the karstic systems. A more noticeable negative trend in spring flow rate is observed for the karst springs characterized by a dominant conduit flow regime and low matrix storage, located in the areas with low cumulative rainfall, and has a stronger relationship with precipitation. Based on the results, decisions on the management of karst water resources should be made considering where the springs bear free surface and pressurized flow conditions.

**Keywords:** climate change; flow condition; new index; spring discharge

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Behrouj Peely, A.; Mohammadi, Z.; Sivelle, V.; Labat, D.; Naderi, M. A New Index to Assess the Effect of Climate Change on Karst Spring Flow Rate. *Sustainability* **2024**, *16*, 1326. https://doi.org/10.3390/su16031326

Academic Editor: Lucio Di Matteo

Received: 23 November 2023

Revised: 7 January 2024

Accepted: 17 January 2024

Published: 4 February 2024

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2024 by the authors.
Licensee MDPI, Basel, Switzerland.
This article is an open access article
distributed under the terms and
conditions of the Creative Commons
Attribution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).

1. Introduction

Climate change, especially water shortage, drought, and floods, might drastically hit both the quantity and quality of aquifers; groundwater has a substantial role in ecosystem preservation and the empowerment of human reconciliation against climate change [1] Karst water resources, because of their vulnerable conduit system or, more generally, because of their fast flow component, are more at risk of facing critical conditions [2]. They are vital groundwater resources worldwide due to the continuously increased water demand for domestic, industrial, and agricultural usage. Based on [3], 9.2% of the world's population depends on drinkable water resources from karst aquifers. More specifically, 40% of the US groundwater resources are provided by karst [4]. In addition, karst aquifer supplies are even more vital for some regions for providing 50% or even more of drinking water such as Austria, the Dinaric region (Europe), Southwest China, Vienna, Damascus, and Rome [5–8].

<footnote n="1">¹ Department of Earth Sciences, Faculty of Sciences, Shiraz University, Shiraz 7146713565, Iran; zmohammadi@shirazu.ac.ir</footnote>

<footnote n="2">² Research Institute for Geo-Hydrological Protection, National Research Council, 06128 Perugia, Italy</footnote>

<footnote n="3">³ Université Toulouse 3-Géosciences Environnement Toulouse-CNRS-UT3-IRD, 31400 Toulouse, France; david.labat@get.omp.eu</footnote>

<footnote n="4">⁴ Department of Earth Sciences, Institute for Advanced Studies in Basic Sciences (IASBS), Zanjan 4513766731, Iran</footnote>

<footnote>* Correspondence: behrouj.a@gmail.com</footnote>
To study the effect of climate change on the karstic springs, it is first necessary to employ a proper approach to identify the projected climate change and to understand climate variability. Downscaling is a technique performed on general circulation models' (GCMs) output to provide climate information at a more suitable scale for local investigation. Empirical relationships between large-scale atmospheric and local climate variables can be expressed with statistical downscaling approaches [9]. Long Ashton Research Station weather generator (LARS-WG) is a stochastic weather generator; it is used to derive statistical variables from observation data in order to create synthetic weather data [10]. This approach has proved efficient for generating climatic variables under climate scenarios, particularly in the calculation of wet/dry spell length [11]. In the most recent report from the Intergovernmental Panel on Climate Change (IPCC), the socio-economic changes and trends in scenarios are described as the shared socio-economic pathway or SSP [12]. There are five illustrative scenarios defined in [12], which cover very low and low emission of greenhouse gases (GHGs) (SSP1-1.9 and SSP1-2.6) to their intermediate (SSP2-4.5), high, and very high emission (SSP3-7.0 and SSP5-8.5). The output of ensemble models of Coupled Model Intercomparison Project 6 (CMIP6) is used for downscaling; CMIP6 explains the earth's response to forcing and an assessment of climate change predictability amongst other descriptions [13].

There are several methods for studying climate change impacts on groundwater flow in karst aquifers including the semi-distributed lumped model [14], integrated hydrological model (MOHISE) [15], soil and water assessment tool (SWAT), and MODFLOW in a one-way connection method [16]. However, these methods need an immense amount of data such as hydraulic conductivity and heterogeneities in physical properties to successfully simulate groundwater flow, which is often extensively challenging in karst aquifers. On the other hand, approaches of lumped parameter modeling on a statistical basis offer solutions with less data (parameter) demand. Ref. [17] used correlation stochastic processes to investigate the rainfall and karst spring relationship. They concluded that the wavelet analyses approach outperforms the linear-based relation of input-output models (i.e., rainfall-runoff) in karst systems which noticeably show non-stationary non-linearity behavior. In other research, continuous and discrete wavelet transforms were used to study the impact of climate change and human activities on karst aquifers [18]. A cross-correlation function (CCF) was used at multi-resolution levels, helping to decompose time series using wavelet analysis. CCF has the potential to be used in karst development analysis. For example, Ref. [19] used cross-correlation with the help of other statistical analyses in time series such as the phase function in two-time series in the Dinaric karst of Croatia and showed that the region has advanced karst. In another study in the Dinaric karst (Rumin Springs), CCF, partial cross-correlation function (PCCF), autocorrelation analysis function (ACF), and spectral density function (SDF) were used in order to study the characteristics of the watersheds adjacent to the karst and their response to rainfall events [20]. Cross-correlation was also used with unprocessed and processed data to identify the contribution of rainfall and time of effect through the main infiltration pathways in karstic aquifers of Campania, southern Italy [21]. In a recent study in southern Italy (Caposele Spring), Ref. [22] introduced a new criterion grounded on the frequency of events surpassing a particular extreme value; they used it to detect acceptable distributions for the data to assess spring discharge extreme events. Time series analysis can also be used to drive rainfall intensity and time lag curves, which reveal valuable information on spring responses and karst segments [23].

Iran, with no exception to the other countries of the world, is experiencing climate change. Ref. [24] estimated higher precipitation in the wet regions (north and west) and lower in the dry regions (south and east), these trends would lead to an increase and decrease in groundwater recharge in these parts of Iran, respectively. Ref. [25] showed that carbonate rocks in Iran cover about 11% of the country and some sites such as Zagros illustrate more development in karstification [26]. To give some sense of the importance of hard rock (i.e., including karst and fractured aquifers) water resources in Iran, it was estimated that 320 and 12,632 m³/year groundwater are estimated to be extracted or
discharged from 920 wells and 37,490 springs (75% of which from carbonate formations), respectively [27]. In addition, hydrochemistry samples from some karst regions in Iran show good quality in terms of low electrical conductivity, often less than 1000 µS/cm [28], which makes them a valuable resource for drinking water supply. So far and to the knowledge of the authors, the impact of climate change on hydrogeological conditions in karstic regions in Iran is rarely reported. An investigation by [29,30] on Bibitarkhon Spring in the Lali region revealed no considerable change in the spring flowrate using an artificial neural network (ANN) with hopeful durability. But in a drier region in south-central Iran (Firouzabad), an evaluation the effect of climate change showed that the quantity and quality of the Firouzabad river and some saline and karstic springs decrease using the coupled global climate model (CGCM) [31].

The characteristics of the karst systems in each part of the Zagros region are highly dependent on the local geological and hydrogeological setting [32]. Therefore, mathematical simulation of the hydrogeological system of its aquifers for spring discharge forecasting can be time consuming and not sufficiently robust to investigate climate conditions out of the range of climate conditions observed within the calibration period. Moreover, Refs. [5,33], among others, already demonstrated that model-based prediction of the effect of climate change on groundwater, particularly on karst aquifers (because of their heterogeneity), experienced uncertainty linked to meteorological variables and land use change. Assuming that the transfer function of karst aquifer is constant over the projected periods (~60 years), statistical approaches can be relevant to investigate the impact of change in input (climate and meteorological forcings) on the output (spring discharge). In this study, the objectives are (a) to use time series analysis and covariance function in the historical period and future to characterize the karst system and (b) to introduce new indices to evaluate the impact of climate change on the karstic springs flow rate.

## 2. Study Area

Iran's plate consists of several structural geology units, as shown in Figure 1 including Zagros, Alborz, Central Iran, and Kopet-Dagh [34]. The Zagros range, located in the Alpine-Himalayan belt, extends over 1600 km from the northwest (in Turkey and Iraq) to the southeast (Oman Sea) in the west of Iran, which contains 40% of Iran's water resources [35]. Fourteen meteorological stations were selected to investigate the climate change impact on karst springs (Figure 2). Masjedsoleiman, with the lowest elevation of 321 m, is the hottest place (mean temperature is 26 °C) among the stations while Kuhrang synoptic station is the coldest place (mean temperature is 10 °C) with the highest elevation located at 2365 m. The minimum and maximum average annual rainfall, with the amount of 279 and 1301 mm/year, occurred at the Fasa and Kuhrang stations in the southern and almost central parts of the region, respectively.

Based on the acquired data from Iran's Ministry of Energy, there are more than 2000 karst springs in eight provinces, whose distribution is shown in Figure 1. The minimum flow of these springs is less than 0.1 to up 1600 L/s and their maximum flow is less than 0.1 to more than 8600 L/s. Moreover, the electrical conductivity varies from 85 to 10,830 µS/cm with an average of 709 µS/cm. The emergence height of these springs varies from 201 to 3850 m above sea level (masl). On the basis of data availability and spatial distribution, 10 springs including Biston, Bernaj, Sarabgarm, Tangsiab, Cureh, Todehzan, Pirgahr, Dimeh, Barm, and Sasan were selected for further study to highlight climate change effects (Table 1 and Figure 2). According to the studies of [36], the reaction time of Sasan spring to maximum rainfall events is from 6 to 85 days with an average of 43 days, its memory effect using ACF is 90 days and, based on CCF, the reaction time for quick flow and base flow is estimated at 18 and 58 days, respectively. It is reported that this spring has a conduit flow type [37]. Rain recording stations were selected for this study based on the closest measuring station of meteorological parameters to the area where each spring appeared (Table 1). In addition, such a station should have enough historical period data
for downscaling. The overall gradient variation in precipitation and flow rate is named as
trend in Table 1.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Karstic springs of Zagros with the selected springs.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** The study area and the distribution of selected stations on a digital elevation map.</figcaption></figure>

To determine the approximate range of the catchment area of each spring, SWAT was used, implementing the digital elevation model (DEM) extracted from the global U.S.
Geological Survey's public domain geographic database HYDRO1k. Finally, taking into consideration other karst springs of each area and also the expansion of karst anticlines, the range from which the springs are fed was determined approximately. The most important geological formations around the springs are as follows: Sasan, Tangsiab, and Sarabgarm Springs: Asmari limestone formation (Miocene); Barm Spring: quaternary sediments and Asmari and Aghajari formations (Miocene), Pirghar and Dime limestones equivalent to the Asmari, Jahorm and Shahbazan formations with Eocene-Oligocene age; Cureh: Permian limestones; Tudezan: Jurassic-Cretaceous limestone and sandstone; and Bernaj and Biston: Biston limestone (Triassic-Cretaceous) and Cretaceous limestones.

<table><caption>**Table 1.** General characteristics of the selected springs with their nearest meteorological station.</caption><thead><tr><th>Spring</th><th>Mean Discharge Q (L/s)</th><th>Catchment Area (km²)</th><th>Q Trend</th><th>Meteorological Station</th><th>Elevation (m)</th><th>Rainfall P (mm/year)</th><th>P Trend</th></tr></thead><tbody><tr><td>Todehzan</td><td>158</td><td>171</td><td>-0.01</td><td>Brojerd</td><td>1629</td><td>473</td><td>0.0003</td></tr><tr><td>Cureh</td><td>218</td><td>99.2</td><td>-0.02</td><td>Brojerd</td><td>1629</td><td>473</td><td>0.0003</td></tr><tr><td>Biston</td><td>674</td><td>34.4</td><td>-0.02</td><td>Polechehr</td><td>1270</td><td>372</td><td>-0.00004</td></tr><tr><td>Tangsiab</td><td>1344</td><td>130</td><td>-0.04</td><td>Tangsiab</td><td>900</td><td>394</td><td>-0.0004</td></tr><tr><td>Sasan</td><td>1686</td><td>351</td><td>-0.05</td><td>Ghaemieh</td><td>922</td><td>548</td><td>-0.0007</td></tr><tr><td>Pirgahr</td><td>1818</td><td>72.1</td><td>-0.12</td><td>Farsan</td><td>2062</td><td>414</td><td>-0.0037</td></tr><tr><td>Sarabgarm</td><td>1824</td><td>62.9</td><td>-0.06</td><td>Sarpolezahab</td><td>545</td><td>422</td><td>-0.0011</td></tr><tr><td>Bernaj</td><td>1869</td><td>193</td><td>-0.03</td><td>Polechehr</td><td>1270</td><td>372</td><td>-0.00004</td></tr><tr><td>Barm</td><td>2183</td><td>556</td><td>-0.05</td><td>Lordgan</td><td>1611</td><td>535</td><td>-0.0009</td></tr><tr><td>Dimeh</td><td>2960</td><td>310</td><td>-0.09</td><td>Kuhrang</td><td>2365</td><td>1309</td><td>-0.0033</td></tr></tbody></table>

### 3. Methodology

The three sequence steps of the study are (1) downscaling the GCMs climate output with LARS-WG6, (2) time series analysis using CCF as well as covariance function, and (3) introducing a new index. The first step is to predict future precipitation; in the second one, the relation between the precipitation and spring flow rate in the historical period is investigated. In the third step, based on such relationship and predicted precipitation, the climate change impact on spring flow rate is studied.

The linear model type was applied to the time series of both target variables (i.e., rainfall and spring discharge) in trend analysis to determine whether such variables have undergone any change. Daily precipitation projection is performed considering two 30-year-long periods (i.e., 2021–2050 for the near future and 2051–2080 for the far future) applying LARS-WG6 with a downscaling technique. The daily output of GCMs in CMIP6 was used under three SSP scenarios including SSP1-1.9, SSP2-4.5, and SSP5-8.5 are representative of low, intermediate, and high emissions of GHGs, respectively. We conducted downscaling for each station with 25 years of observed data from 1990 to 2014. The correlation coefficient between the observed and generated data for all stations was above 0.98 except for Sarpolezahab, which was 0.96.

The CCF can be used to determine the degree of similarity between two signals or two components. The normalized CCF $r_{xy}(k)$ of two time series $x$ and $y$ is calculated as follows:

$$r_{xy}(k) = \begin{cases} \frac{C_{xy}(k)}{\sigma_x \sigma_y}; & k \ge 0 \\ \frac{C_{yx}(-k)}{\sigma_x \sigma_y}; & k < 0 \end{cases} \qquad (1)$$

where $k$ is the time lag, $\sigma_x$ and $\sigma_y$ are the standard deviation of $x$ or rainfall and $y$ or spring flow rate, respectively. When $n$ is the length of the time series and $\mu_x$ and $\mu_y$ are the average of $x$ and $y$, respectively, $C_{xy}(k)$ is calculated from the following equation:

$$C_{xy}(k) = \frac{1}{n} \sum_{t=1}^{n-k} (x_t - \mu_x)(y_{(t+k)} - \mu_y) \qquad (2)$$
The correlation between the elements of a series with other elements of the same series ($r_{yy}(k)$) that are separated by a certain time interval is widely used in ACF in time series investigations of karst springs:

$$r_{yy}(k) = \frac{C_{yy}(k)}{\sigma_y^2}; k \ge 0 \qquad (3)$$

The covariance function ($C_{xy}(k)$) explains how much two random variables with different spatial or temporal separation change and describe the second-order dependence of random processes. The popularity of covariance functions in spatial and space-time statistics is due to the fact that the properties of Gaussian random fields are fully described by first- and second-order moments. Therefore, covariance functions are very important for the estimation and prediction of Gaussian random fields. This assumes that the value of the time delay calculated in the historical period will be the same in the future, which is a correct assumption considering that the karst system in terms of conduit characteristics will not change during the next 60 years. Therefore, $C_{xy}(k)$ or $C_{PQ}(k)$, which is obtained based on the measured data and is a function of the covariance between rainfall and the springs flow rate, does not change in the future. Hence, the standard deviation of the spring flow rate ($\sigma_Q$) in the future can be calculated from the following formula:

$$\sigma_Q = \frac{C_{PQ}(k)}{\sigma_P r_{PQ}(k)} \qquad (4)$$

where $\sigma_P$ is the standard deviation of the rainfall and can be calculated from the predicted rainfall data (simulated based on climate models and downscaling at the site of a meteorological station) and $r_{PQ}(k)$ is the correlation coefficient between rainfall and the springs flow rate, which is obtained using the data of the historical period. Considering $N$ as the number of measurements, $Q_t$ is the spring flow rate at any time and $Q_{mean}$ is the average flow rate of the springs:

$$\sigma_Q = \sqrt{\frac{\sum_{t=1}^{N} (Q_t - Q_{mean})^2}{N}} \qquad (5)$$

Therefore, it is possible to calculate the sum of the deviation of the flow rate of springs from their average value based on the following formula:

$$\sum_{t=1}^{N} (Q_t - Q_{mean})^2 = N \left( \frac{C_{PQ}(k)}{\sigma_P r_{PQ}(k)} \right)^2 \qquad (6)$$

Then, by using the following two formulas, the effect of climate change on the flow rate of the springs can be assessed as follows (in fact, the sum of the deviation of the flow rate of the springs from their average value in the historical period (b) and the future (f) are compared):

$$dQ_d = \sum_{t=1}^{N} (Q_t - Q_{mean})_f^2 - \sum_{t=1}^{N} (Q_t - Q_{mean})_b^2 \qquad (7)$$

For dimensionless formulas, the effect can be calculated as follows:

$$I_{Q_d} = \frac{\sum_{t=1}^{N} (Q_t - Q_{mean})_f^2}{\sum_{t=1}^{N} (Q_t - Q_{mean})_b^2} \qquad (8)$$

The following criterion shows differentiation of the changes in the flow rate of the springs from the average value in the historical period and future periods compared to the historical period:

$$I_{dQ_d} = \frac{dQ_d}{\sum_{t=1}^{N} (Q_t - Q_{mean})_b^2} \qquad (9)$$

Since $\sum_{t=1}^{N} (Q_t - Q_{mean})^2$ is too large and for the ease of comparison, it is square rooted and, hereinafter, is referred to as $d$. Finally, the following formula is introduced as a new
criterion for evaluating the effect of climate change on the karstic springs discharge. In this criterion, the differentiation of the changes in the flow rate of the springs from the average value in the historical period and future periods are representative of the hydrogeological system (karst) and the precipitation changes in the same two periods are representative of the precipitation variability under the influence of climate change:

$$I_{cc} = (|I_{dQ_d}|) \left(\frac{dP}{P_b}\right) \qquad (10)$$

$dP = P_f - P_b$ is the difference between the precipitation in the future climate change scenario ($P_f$) and the historical period ($P_b$).

Table 2 shows the values of $I_{Q_d}$, $I_{dQ_d}$ and $I_{cc}$ with the corresponding possible conditions of the groundwater flow in karst. $I_{cc}$ is somewhat different from the other criteria; therefore, the corresponding values and springs covered in such values are also different. The terms including “overflow, flooding, pressurized and free surface flow” in Table 2 are based on those in the literature, such as in the work of [38–40], and adapted according to their use in this study. The value of $I_{Q_d}$ is less than 0.50 (when $\sum_{t=1}^{N}(Q_t - Q_{mean})_b^2$ is more than twice $\sum_{t=1}^{N}(Q_t - Q_{mean})_f^2$) and also, the two criteria, $I_{dQ_d}$ less than 0.50 and $I_{cc}$ more than 0.25, show that the underground water flow is beyond the capacity of the conduits and the higher water head in the conduits in conditions of high flow gradient cause the matrix to be fed (flooding/overflow). In addition, the excess amount of rainwater flows in the form of runoff outside the capacity of the entrance of the conduits. $I_{cc}$ in the range of 0.10–0.25 (equivalent to $I_{dQ_d}$ [−0.50, −0.20] and $I_{Q_d}$ 0.50–0.80) can indicate that the underground water flows in the entire section of the conduits and can still charge the matrix (pressurized flow conditions). $I_{cc}$ in the interval 0.05–0.10 shows the amount of groundwater flow increases mildly. In the two intervals of 0.00–0.05 and (−0.05)–0.00, there is no significant change in the underground water flow rate and the conditions are almost similar to those before the climate change scenario; in these two intervals, the flow rate of the springs increases and decreases slightly, respectively. $I_{cc}$ between −0.10 and −0.05 shows groundwater flow decreases mildly. $I_{cc}$ in the range of (−0.25)–(−0.10) can indicate that the groundwater flows only in a part of the conduit section (free surface flow conditions) and the conduit is fed from the matrix. When $I_{Q_d}$ is greater than 1.50 ($\sum_{t=1}^{N}(Q_t - Q_{mean})_f^2$ is more than one and a half times $\sum_{t=1}^{N}(Q_t - Q_{mean})_b^2$), $I_{dQ_d}$ greater than 0.50 and $I_{cc} < (−0.25)$ indicates that the groundwater flow is less than the capacity of the conduits and the extremely low head of water in the conduits in the condition of very low flow gradient causes matrix discharge. In this situation, if the spring is located in a low rainfall area and its dominant flow regime is a conduit type, there is a high possibility of drying up. Figure 3 demonstrates a conceptual idealized model of climate change conditions in a karst system with a doline inlet and a spring outlet with typical formations in Zagros including Gachsaran evaporite (Gs), Asmari karstic limestone (As), and impermeable Pabdeh (Pb) and Gurpi (Gp). The lithological units are shale and marl Pabdeh–Gurpi formations (Pb-Gu: Santonin-Oligocene), Asmari-Jahrom limestone formations (As-Ja: Pliocene–Miocene), Gachsaran gypsum formation (Gs: Oligocene–Miocene), and lakes and rivers deposits (Q) of present-day. More information about these formations can be found in [36] studying Sasan spring. In Figure 3, the conditions of the water table, conduits, and karstic spring in four main states under the influence of climate change are illustrated: (a) dry condition of the conduit and spring, (b) conduit with free-surface water, (c) conduit full of water and under pressure, and (d) flooding conditions. The other four situations discussed in Table 2 are classified under conditions b and c depending on the climatic changes. If there is a slight decrease or increase in the flow, the conditions of the karst aquifer will approach b and c states, respectively. Two states, a and d, occur in extreme events: (a) in the case of a sharp decrease in precipitation and (d) in the case of a sharp increase in precipitation.
<table><caption>**Table 2.** The range of *I*<sub>Qd</sub>, *I*<sub>dQd</sub>, and *I*<sub>cc</sub> and the concept of possible groundwater flow conditions in idealized karst aquifers are illustrated in Figure 3.</caption><thead><tr><th>I<sub>Qd</sub></th><th>I<sub>dQd</sub></th><th>I<sub>cc</sub></th><th>Possible Corresponding Flow Conditions in Karst</th><th>Description</th></tr></thead><tbody><tr><td>&lt;0.50</td><td>&lt;(-0.50)</td><td>&gt;0.25</td><td>Overflow/Flooding</td><td>Groundwater flow is beyond conduit capacity; it recharges the matrix and the excess amount flows as runoff (back-flooding).</td></tr><tr><td>0.50-0.80</td><td>(-0.50)-(-0.20)</td><td>0.10-0.25</td><td>Pressurized flow</td><td>Groundwater flows in the whole cross section of conduits and it recharges the matrix.</td></tr><tr><td>0.80-0.90</td><td>(-0.20)-(-0.10)</td><td>0.05-0.10</td><td>Mild flow increase</td><td>The growth of groundwater flow is mild.</td></tr><tr><td>0.90-1</td><td>(-0.10)-0.00</td><td>0.00-0.05</td><td>Little flow increase</td><td>There is no much difference with the previous flow conditions apart from a trivial flow rise.</td></tr><tr><td>1-1.10</td><td>0.00-0.10</td><td>(-0.05)-0.00</td><td>Little flow decline</td><td>There is no much difference with the previous flow conditions apart from a trivial flow reduction.</td></tr><tr><td>1.10-1.20</td><td>0.10-0.20</td><td>(-0.05)-(-0.10)</td><td>Mild flow decline</td><td>The reduction in groundwater flow is mild</td></tr><tr><td>1.20-1.50</td><td>0.20-0.50</td><td>(-0.10)-(-0.25)</td><td>Free surface flow</td><td>Groundwater flows in some part of conduits and it discharges the matrix.</td></tr><tr><td>&gt;1.50</td><td>&gt;0.50</td><td>&lt;(-0.25)</td><td>Spring dryness</td><td>The possibility of spring going dry is highly likely, especially in the event of low average rainfall and dominancy of conduit flow; soil-matrix flow may become dominant in case of soil layer existence.</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Conceptual model of climate change impacts of hydrological functioning of a karst system: (a) dry condition of conduit and spring, (b) conduit with free water surface, (c) flooded conduit, and (d) flood condition. New Age sediments (Q), Gachsaran formation (Gs), Asmari (As), Pabdeh (Pb) and Gurpi (Gp).</figcaption></figure>

# 4. Results and Discussion

## 4.1. Precipitation under Climate Change

Precipitation increases in the two 30-year time periods of the future for some places even to an amount of 340 mm with $dP/P_b$ of about 0.55 under the SSP1-1.9 scenario
(Figure 4). The amount of precipitation multiplies by about 0.50 in Tangsiab in the first time period and Farsan in the second time period. The worst possible scenario in terms of precipitation reduction for the Zagros region in this study is SSP2-4.5. At stations Gaemieh and Kuhrang, $dP$ is in the range of -80 to -120 mm/year and -240 to -280 mm/year with $dP/P_b$ of -0.15 to -0.22, respectively. Tangsiab, Polechehr, and Brojerd experienced almost the same rainfall reduction ($dP = 40$ mm/year) (Figure 4). Likewise, Farsan and Lordgan do not tolerate much precipitation reduction in the near future (almost 18 to 27 mm/year) but in the far future, $dP$ is about 55–67 mm/year with $dP/P_b$ of -0.09 to -0.12. In the first 30-year period, future precipitation decreases to a maximum of 22 mm and the ratio of $dP/P_b$ is about -0.04 at Sarpolezahab and Polechehr stations under the SSP5-8.5 scenario. However, $dP/P_b$ decreases more to about -0.08 for some places such as Kuhrang and Lordgan in the second 30-year period of the future under this scenario. Brojerd and Tangsiab experience more precipitation of up to 69 and 50 mm/year than $P_b$ in the near and far periods, respectively.

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Precipitation change ($dP$) with respect to the historical period precipitation ($P_b$) under three scenarios in the two 30-year periods from 2021 to 2080.</figcaption></figure>

Generally, precipitation increases under the SSP1-1.9 scenario and decreases under SSP2-4.5 while it both decreases and increases under the SSP5-8.5 scenario. The northern part, in terms of precipitation reduction, is affected by climate change to a lesser extent than the southern part under SSP2-4.5. In contrast, under SSP5-8.5, the southern part of Zagros suffers no precipitation reduction. Kuhrang, in the central parts, is affected the most, particularly from 2051–2080 under SSP2-4.5. On the whole, the annual average precipitation of all 14 stations in the baseline is 540 mm/year but in the first and second 30-year periods under SSP1-1.9, it increases to 665 and 691 mm/year; under SSP2-4.5, it endures reduction and becomes 504 and 484 mm/year and under SSP5-8.5, it experiences a 13 mm/year rise (553 mm/year) and a 6 mm/year drop (534 mm/year), respectively. The corresponding percentages in the first and second 30-year periods under SSP1-1.9 are 24% and 28% increases, under SSP2-4.5 they are 6% and 9% reductions, and under SSP5-8.5 they are a 3% increase and less than 1% reduction, respectively.

## 4.2. *Time Series Analysis*

The Karst system memory effect is estimated as the time lag where the ACF value first appears lower than 0.2 [41]. Karst systems with a longer memory effect are known as low
karstified or with a large storage volume. A developed karst system, with well-established Karst features such as conduits, is generally characterized by a low memory effect. As observed in Figure 5, the memory effect in Cureh is the highest value of 42 months and, after that, are Barm, Tudehzan, and Sarabgarm Springs with 33.2, 27.5, and 21 months, respectively. The memory effect in the rest of the springs is less than 6 months; among these springs, Sasan has the greatest memory effect at 5.2 months and Pirghar at 1.9 months has the least. It is expected that springs with higher karstification and more conduit regime systems show more correlation with rainfall. From this point of view, as shown in Figure 6, the Biston, Bernaj, and Pirghar Springs have the highest coefficient with values between 0.45 and 0.49. However, Tudehzan, Barm, Cureh, and Sarabgarm Springs have a coefficient of less than 0.25 (from 0.17 to 0.23).

<figure><img src="image_8.png" /><figcaption>**Figure 5.** Chart of ACF values in springs based on monthly temporal resolution: (a) Sasan; (b) Pirghar, Barm, and Dimeh; (c) Tangsiab, Todehzan, and Cureh; and (d) Biston, Bernaj, and Sarabgarm.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Chart of CCF values in springs. (a) Sasan Spring and Gaemieh station. (b) Pirghar Spring and Farsan station, Barm Spring and Lordegan station, and Dimeh Spring and Kohrang station. (c) Tangsiab Spring and Tangsiab station, Tudehzan and Cureh Springs and Borujerd station. (d) Biston and Bernaj Springs and Sararoud station and Sarabgarm Spring and Sarpolezahab station.</figcaption></figure>
4.3. Flow Rate under Climate Change

Figure 7 shows the sum of differences between spring discharge at any time and its average (d) in the historical period (baseline), in the near future (2050–2021), and in the far future (2080–2051) under SSP1-1.9, SSP2-4.5, and SSP5-8.5, respectively. Considering the values of d observed in the historical period, the springs can be divided into three groups: Todehzan and Cureh (changes from 1164 up to 1969 L/s); Sasan, Sarabgarm, Dimeh, Barm, and Tangsiab (changes from 5580 to up 8657 L/s); and the third group Biston, Pirghar, and Bernaj (changes from 11,425 to up 25,715 L/s). The higher value of d means that the karstic system is more active in the water transmission or is more conduit type. This derivation is partially consistent with the trend analysis of changes in the springs’ flow rate; for instance, Pirghar, with a decrease of –0.085, has the largest decrease in flow rate and Tudezan, with a decrease of –0.018, has experienced the lowest decrease in the flow rate. Therefore, the two springs of Pirghar and Bernaj, which show the most changes in the historical period and in all the future scenarios, are the most karstified springs contrasted to the two springs of Tudezan and Cureh, which have the lowest degree of karstification. The other springs are in the range between these two end-members of a low and high degree of karstification (low to high karstified springs: Sasan, Sarabgarm, Dimeh, Barm, Tangsiab, and Biston). Moreover, this conclusion is confirmed by the memory effect and CCF, as discussed in the previous section. The Biston, Bernaj, and Pirghar Springs have the highest CCF coefficient and the lowest memory effect, while this is the opposite in springs like Todehzan and Cureh. The value of d in Bernaj Spring is 25,715 L/s in the historical period; its lowest (23,843 L/s) and highest (38,946 L/s) values occur in the far future under SSP1-1.9 and SSP2-4.5. The changes in Dimeh and Sasan Springs under the influence of climate change scenarios are significant. Dimeh is the only spring whose value of d is higher than the historical period (7880 L/s) in all climate change scenarios.

<figure><img src="image_10.png" /><figcaption>Figure 7. Changes in *d* in springs in the historical period and under three climate change scenarios SSP1-1.9, SSP2-4.5 and SSP5-8.5 in two time periods of the near future (2050–2021) and far future (2051–2080).</figcaption></figure>

Figure 8 shows $dQ_d$ in springs under SSP1-1.9, SSP2-4.5 and SSP5-8.5 in near future and far future. The negative sign indicates that the average discharge of the spring has increased under the influence of the climate change scenario in the future, which happened only in the SSP1-1.9 scenario. Under this scenario, Tudehzan and Cureh Springs have the least changes in both 30-year periods with $dQ_d$ of about (–85 L/s) and (–140 L/s),
respectively. Under the SSP2-4.5 climate change scenario, Barnaj experiences the most changes in near future (8575 L/s) and far future (13,231 L/s). During these two periods, $dQ_d$ is the lowest in Tudehzan Spring and has values of 440 and 544 L/s, respectively. The most $dQ_d$ changes under the SSP5-8.5 scenario occur in the near future in Pirghar Spring with 8812 L/s and in the far future in Bernaj Spring with 10,023 L/s. Like other scenarios, Tudehzan Spring has the least changes, so that $dQ_d$ is about 300 L/s in both periods. Based on $dQ_d$, the Tudehzan, Cureh, Barm and Sarabgarm Springs undergo the least changes and Dimeh, Pirghar and Bernaj bear the most changes. The results of this index show that the effect of the karst system on the rate of changes in the flow of springs is more specific than the changes in precipitation.

Figure 9 shows the changes in $I_{dQ_d}$ in springs under three climate change scenarios, SSP1-1.9, SSP2-4.5, and SSP5-8.5, in the near and far future. The negative sign of this index indicates that the average discharge of the spring has increased under the influence of the climate change scenario in the future, which only happened in the SSP1-1.9 scenario in both time periods. In the near future of the SSP1-1.9 scenario, the average flow rate has decreased by 0.0003, 0.02, 0.11, and 0.16 in the four springs of Biston, Bernaj, Pirghar, and Dimeh, respectively. In the far future, such a condition prevails only for Dimeh Spring with a value of 0.097. Under the SSP2-4.5 scenario, the largest changes in $I_{dQ_d}$ are observed, so that its highest values in both 30-year periods occurs in Dimeh Spring, with values of 0.70 and 0.81. The spring that experiences the least change in both periods is Barm, whose $I_{dQ_d}$ value is 0.17 and 0.21 in the near and far periods, respectively. Under the SSP5-8.5 scenario, in the near and far future, Pirghar Spring tolerates the greatest change and the $I_{dQ_d}$ value is 0.51 and 0.57, respectively. The lowest value of $I_{dQ_d}$ happens in Tangsiab (0.16) and Sarabgarm (0.05) Springs in the near and far future, respectively. In fact, by using this index, it is possible to identify extreme events in the value of the flow rate. The more positive value of this index indicates that there have been more discharges lower than the average value and vice versa. In conclusion, based on this index, most of the springs undergo free surface conditions in the future. Moreover, Dimeh, Sasan, Pirgahr, and Bernaj might even bear dryness (Figure 9).

$I_{cc}$ changes are shown in Figure 10 in springs under SSP1-1.9, SSP2-4.5, and SSP5-8.5 in two time periods of the near future and far future. In this index, in order to investigate the effect of climate change, the direct effect of precipitation changes is considered (Equation (10)). The negative sign of this index indicates that the precipitation has decreased under the influence of the climate change scenario in the future. Under the SSP1-1.9 scenario, Tangsiab Spring shows the greatest increase in changes in both 30-year periods. In this condition, the conduit of this spring is completely filled and it feeds the matrix. Although it has the least changes in the near future, Biston Spring has significant incremental changes of 0.05 in the far future. Sarabgarm and Sasan Springs are other springs that have significant incremental changes in the first 30 years, which is more due to changes in precipitation than the karst system of the spring. Barm Spring shows little changes in both periods under SSP1-1.9. Under the SSP2-4.5 scenario, Dimeh and Sasan Springs bear the most decreasing changes, so that $I_{cc}$ decreases by 0.09 and 0.06 in the near future and by 0.17 and 0.12 in the far future, respectively. Under such circumstances, the free-surface flow is in place and the transfer flow is from the matrix into the conduit. Due to the relatively high average rainfall and the matrix reserve (big catchment area), the possibility of drying up of these springs is low; however, their flow rate decreases considerably. But the Barm Spring in the first 30 years with a decrease in 0.0008 and the Sarabgarm Spring with a decrease in 0.02 in the second 30 years have the least decreasing changes under SSP2-4.5. In the near future, Pirghar Spring has the most increasing changes with 0.027 and Sarabgarm Spring has the most decreasing changes with 0.013 under the SSP5-8.5 scenario. In the far future, Dimeh and Bernaj Springs bear the most decreasing changes with $I_{cc}$ of 0.043 and 0.018, respectively. In closing, based on this index, no springs bear the extreme events of overflooding or drying up; rather, they experience little to mild changes.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** *dQ*<sub>d</sub> value changes in terms of L/s in springs based on climate change scenarios SSP1-1.9, SSP2-4.5 and SSP5-8.5 in two time periods of the near future (2050–2021) and distant (2080–2051).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** *I*<sub>dQ<sub>d</sub></sub> in the springs under SSP1-1.9, SSP2-4.5, and SSP5-8.5 in two time periods of the near future (2050–2021) and far future (2080–2051). See Figure 3 for the concept of different terms of flow regime.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 10.** *I*<sub>cc</sub> in springs under SSP1-1.9, SSP2-4.5, and SSP5-8.5 in two time periods of the near future (2050–2021) and far future (2080–2051). See Figure 3 for the concept of different terms of flow regime.</figcaption></figure>

## 4.4. Limitations and Uncertainties

Data requirement in long term measurements and the interval of monitoring particularly for karst spring discharge is a challenge. In fact, in order to calculate the covariance function so that it is safe to assume that it stays constant in future, such measures are required. On the other hand, the assumption on which the covariance function between rainfall and karst system stays constant in the near future is conceivable but it is uncertain in far future. Another constrain is the necessity that the weather station site should be within the spring catchment area. Moreover, the downscaling technique and the number of global circulation models (GCMs) are other uncertainty sources of such studies. Although the multi-model ensemble of coupled model inter-comparison project phase 6 (CMIP6) were used to resolve the shortcomings of an individual GCM, the adequate number of GCMs is unknown; see [42,43] for more details. To overcome such limitations and uncertainties, the springs and stations with enough data having acceptable monitoring intervals and 37 GCMs were selected. Quantification of the extent of such uncertainties with mitigation strategies is suggested for future research.

# 5. Conclusions

In this study, three indices (i.e., *dQ*<sub>d</sub>, *I*<sub>dQd</sub>, and *I*<sub>cc</sub>), taking the covariance function properties into account, were introduced to investigate the impact of climate change on karstic spring discharge. *dQ*<sub>d</sub> shows the variability of spring discharge from past to future; *I*<sub>dQd</sub> uses such variability over the historical data and *I*<sub>cc</sub> considers the effect of precipitation and spring discharge change together. The feasible groundwater flow
conditions were classified between two extremes, namely from $I_{cc} > 0.25$ to $I_{cc} < (-0.25)$ regarding the indication of overflow/flooding and spring dryness; pressurized and free surface flow is between such extremes. The results were assessed and validated with the karstic characteristics of the springs so that the springs with high matrix storage and less conduit flow such as Todehzan and Cureh are less impacted by climate change; whereas, the springs with dominant conduit flow regime such as Pirgahr and Bernaj suffer more. In addition, there are some other types of springs such as Dimeh and Sasan that, because of significant precipitation change, despite having a great catchment area, likely bear free surface conditions. Overall, except for Dimeh and Sasan under the worst climate change scenario (SSP2-4.5), which might undergo dryness, all the springs experience the flow conditions between pressurized and free surface flow. The introduced indices can be used for any karstic spring in the world to time-economically assess the impact of climate change on the flow rate for water resources management.

**Author Contributions:** Conceptualization, A.B.P. and Z.M.; Software, A.B.P. and M.N.; Validation, Z.M., V.S. and D.L.; Formal analysis, A.B.P.; Investigation, A.B.P. and M.N.; Data curation, A.B.P. and M.N.; Writing—original draft, A.B.P.; Writing—review & editing, Z.M., V.S. and D.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Iran National Science Foundation (INSF), under Grant Number 99011851.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

**Abbreviations**

<table><tr><td>ACF</td><td>autocorrelation analysis function</td></tr><tr><td>ANN</td><td>artificial neural network</td></tr><tr><td>CCF</td><td>cross correlation function</td></tr><tr><td>CGCM</td><td>coupled global climate model</td></tr><tr><td>CMIP6</td><td>coupled model intercomparison project 6</td></tr><tr><td>DEM</td><td>digital elevation model</td></tr><tr><td>GCM</td><td>general circulation models</td></tr><tr><td>GHGs</td><td>greenhouse gases</td></tr><tr><td>IPCC</td><td>intergovernmental panel on climate change</td></tr><tr><td>LARS-WG</td><td>Long Ashton Research Station weather generator</td></tr><tr><td>masl</td><td>meters above sea level</td></tr><tr><td>PCCF</td><td>partial cross correlation function</td></tr><tr><td>SDF</td><td>spectral density function</td></tr><tr><td>SSP</td><td>shared socio-economic pathway</td></tr><tr><td>SWAT</td><td>soil and water assessment tool</td></tr><tr><td>Notations</td><td></td></tr><tr><td>C<sub>PQ</sub>(k)</td><td>the covariance function between rainfall and the springs flow rate</td></tr><tr><td>C<sub>xy</sub>(k)</td><td>the covariance function</td></tr><tr><td>d</td><td>square root of Σ<sup>N</sup><sub>t=1</sub>(Q<sub>t</sub> - Q<sub>mean</sub>)<sup>2</sup></td></tr><tr><td>dP</td><td>is the difference between the precipitation in the future climate change scenario (P<sub>f</sub>) and the historical period (P<sub>b</sub>)</td></tr><tr><td>dQ<sub>d</sub></td><td>variability of spring discharge from past to future</td></tr><tr><td>I<sub>Qd</sub></td><td>variability index of spring discharge from past to future</td></tr><tr><td>I<sub>dQd</sub></td><td>spring discharge variability over the historical data</td></tr><tr><td>I<sub>cc</sub></td><td>effect of precipitation and spring discharge change together</td></tr><tr><td>r<sub>PQ</sub>(k)</td><td>the correlation coefficient between rainfall and the springs flow rate</td></tr></table>
<table><tr><td>r<sub>yy</sub>(k)</td><td>correlation between the elements of a series with other elements of the same series</td></tr><tr><td>k</td><td>time lag</td></tr><tr><td>n</td><td>the length of a time series</td></tr><tr><td>N</td><td>the number of measurements</td></tr><tr><td>P</td><td>rainfall</td></tr><tr><td>P<sub>b</sub></td><td>precipitation in the historical period</td></tr><tr><td>P<sub>f</sub></td><td>precipitation in the future climate change scenario</td></tr><tr><td>Q</td><td>spring flow rate</td></tr><tr><td>Q<sub>mean</sub></td><td>the average flow rate of the springs</td></tr><tr><td>Q<sub>t</sub></td><td>the spring flow rate at any time</td></tr><tr><td>μ<sub>x</sub></td><td>the average of x</td></tr><tr><td>μ<sub>y</sub></td><td>the average of y</td></tr><tr><td>σ<sub>x</sub></td><td>the standard deviation of x</td></tr><tr><td>σ<sub>y</sub></td><td>the standard deviation of y</td></tr><tr><td>σ<sub>p</sub></td><td>the standard deviation of the rainfall</td></tr><tr><td>σ<sub>Q</sub></td><td>the standard deviation of the springs flow rate</td></tr></table>

## References

1. Taylor, R.G.; Scanlon, B.; Döll, P.; Rodell, M.; Van Beek, R.; Wada, Y.; Longuevergne, L.; Leblanc, M.; Famiglietti, J.S.; Edmunds, M. Ground water and climate change. *Nat. Clim. Chang.* **2013**, *3*, 322–329. [CrossRef]

2. Hartmann, A.; Jasechko, S.; Gleeson, T.; Wada, Y.; Andreo, B.; Barberá, J.A.; Brielmann, H.; Bouchaou, L.; Charlier, J.B.; Darling, W.G.; et al. Risk of groundwater contamination widely underestimated because of fast flow into aquifers. *Proc. Natl. Acad. Sci. USA* **2021**, *118*, e2024492118. [CrossRef] [PubMed]

3. Stevanović, Z. Karst waters in potable water supply: A global scale overview. *Environ. Earth Sci.* **2019**, *78*, 662. [CrossRef]

4. Quinlan, J.F.; Ewers, R.O. Subsurface drainage in the Mammoth Cave area. In *Karst Hydrology*, Springer: Boston, MA, USA, 1989; pp. 65–103.

5. Hartmann, A.; Goldscheider, N.; Wagener, T.; Lange, J.; Weiler, M. Karst water resources in a changing world: Review of hydrological modeling approaches. *Rev. Geophys.* **2014**, *52*, 218–242. [CrossRef]

6. Al-Charideh, A. Recharge rate estimation in the Mountain karst aquifer system of Fige h spring, Syria. *Environ. Earth Sci.* **2012**, *65*, 1169–1178. [CrossRef]

7. Kresic, N.; Stevanovic, Z. (Eds.) *Groundwater Hydrology of Springs: Engineering, Theory, Management and Sustainability*; Butterworth-Heinemann: Waltham, MA, USA, 2009.

8. Wu, P.; Tang, C.; Zhu, L.; Liu, C.; Cha, X.; Tao, X. Hydrogeochemical characteristics of surface water and groundwater in the karst basin, southwest China. *Hydrol. Process.* **Int. J. 2009**, *23*, 2012–2022. [CrossRef]

9. Gutiérrez, J.M.; San-Martín, D.; Brands, S.; Manzanas, R.; Herrera, S. Reassessing statistical downscaling techniques for their robust application under climate change conditions. *J. Clim.* **2013**, *26*, 171–188. [CrossRef]

10. Semenov, M.A.; Barrow, E.M. *A Stochastic Weather Generator for Use in Climate Impact Studies*; User Man Herts: Hertfordshire, UK, 2002.

11. Hassan, Z.; Shamsudin, S.; Harun, S. Application of SDSM and LARS-WG for simulating and downscaling of rainfall and temperature. *Theor. Appl. Climatol.* **2014**, *116*, 243–257. [CrossRef]

12. IPCC. Summary for Policymakers. In *Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change*; Masson-Delmotte, V., Zhai, P., Pirani, A., Connors, S.L., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M.I., et al., Eds.; Cambridge University Press: Cambridge, UK, 2021.

13. Eyring, V.; Bony, S.; Eyring Meehl, G.A.; Senior, C.A.; Stevens, B.; Stouffer, R.J.; Taylor, K.E. Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization. *Geosci. Model Dev.* **2016**, *9*, 1937–1958. [CrossRef]

14. Dubois, E.; Doummar, J.; Pistre, S.; Larocque, M. Calibration of a lumped karst system model and application to the Qachouch karst spring (Lebanon) under climate change conditions. *Hydrol. Earth Syst. Sci.* **2020**, *24*, 4275–4290. [CrossRef]

15. Brouyère, S.; Carabin, G.; Dassargues, A. Climate change impacts on groundwater resources: Modelled deficits in a chalky aquifer, Geer basin, Belgium. *Hydrogeol. J.* **2004**, *12*, 123–134. [CrossRef]

16. Pisinaras, V. Assessment of future climate change impacts in a Mediterranean aquifer. *Glob. NEST J.* **2016**, *18*, 119–130.

17. Labat, D.; Ababou, R.; Mangin, A. Rainfall–runoff relations for karstic springs. Part II: Continuous wavelet and discrete orthogonal multiresolution analyses. *J. Hydrol.* **2000**, *238*, 149–178. [CrossRef]

18. Charlier, J.B.; Ladouche, B.; Maréchal, J.C. Identifying the impact of climate and anthropic pressures on karst aquifers using wavelet analysis. *J. Hydrol.* **2004**, *523*, 610–623. [CrossRef]

19. Pavlić, K.; Parlov, J. Cross-correlation and cross-spectral analysis of the hydrographs in the northern part of the Dinaric karst of Croatia. *Geosciences* **2019**, *9*, 86. [CrossRef]

20. Denić-Jukić, V.; Lozić, A.; Jukić, D. An Application of Correlation and Spectral Analysis in Hydrological Study of Neighboring Karst Springs. *Water* **2020**, *12*, 3570. [CrossRef]
21. Fiorillo, F.; Doglioni, A. The relation between karst spring discharge and rainfall by cross-correlation analysis (Campania, southern Italy). *Hydrogeol. J.* **2010**, *18*, 1881–1895. [CrossRef]

22. Leone, G.; Pagnozzi, M.; Catani, V.; Ventafridda, G.; Esposito, L.; Fiorillo, F. A hundred years of Caposele spring discharge measurements: Trends and statistics for understanding water resource availability under climate change. *Stoch. Environ. Res. Risk Assess.* **2021**, *35*, 345–370. [CrossRef]

23. Martín-Rodríguez, J.F.; Mudarra, M.; De la Torre, B.; Andreo, B. Towards a better understanding of time-lags in karst aquifers by combining hydrological analysis tools and dye tracer tests. Application to a binary karst aquifer in southern Spain. *J. Hydrol.* **2023**, *621*, 129643. [CrossRef]

24. Abbaspour, K.C.; Faramarzi, M.; Ghasemi, S.S.; Yang, H. Assessing the impact of climate change on water resources in Iran. *Water Resour. Res.* **2009**, *45*, 1157–1172. [CrossRef]

25. Raeisi, E.; Kowsar, N. Development of Shahpour Cave, southern Iran. *Cave Karst Sci.* **1997**, *24*, 27–34.

26. Vardanjani, H.K.; Bahadorinia, S.; Ford, D.C. An Introduction to Hypogene Karst Regions and Caves of Iran. In *Hypogene Karst Regions and Caves of the World*; Springer: Cham, Switzerland, 2017; pp. 479–494.

27. Rostam Afshar, N.; Kazemi, H.; Nobahar, F. Quality Protection Legislations for Karst Water Resources. *Iran-Water Resour. Res.* **2010**, *5*, 56–58.

28. Minoyi, M. The Role of Fractures in the Flow of Underground Water in the Karst Region of Shahu Mountains, Kurdistan. Master's Thesis, Shahrood University of Technology, Shahrud, Iran, 2010.

29. Zeydalinejad, N.; Nassery, H.R.; Alijani, F.; Shakiba, A. Forecasting the resilience of Bibitarkhoun karst spring, southwest Iran, to the future climate change. *Model. Earth Syst. Environ.* **2020**, *6*, 2359–2375. [CrossRef]

30. Zeydalinejad, N.; Nassery, H.R.; Shakiba, A.; Alijani, F. Prediction of the karstic spring flow rates under climate change by climatic variables based on the artificial neural network: A case study of Iran. *Environ. Monit. Assess.* **2020**, *192*, 375. [CrossRef]

31. Naderi, M.; Raeisi, E.; Zarei, M. The impact of halite dissolution of salt diapirs on surface and ground water under climate change, South-Central Iran. *Environ. Earth Sci.* **2016**, *75*, 708. [CrossRef]

32. Peely, A.B.; Mohammadi, Z.; Raeisi, E. Breakthrough curves of dye tracing tests in karst aquifers: Review of effective parameters based on synthetic modeling and field data. *J. Hydrol.* **2021**, *602*, 126604. [CrossRef]

33. Riedel, T.; Weber, T.K. The influence of global change on Europe's water cycle and groundwater recharge. *Hydrogeol. J.* **2020**, *28*, 1939–1959. [CrossRef]

34. Aghanabati, A. *Geology of Iran*; Organization of Geology and Mineral Explorations of Iran: Tehran, Iran, 2004.

35. Moloudia, F.; Shokati, S. *Assessment of Climate Change in Northern Zagros Forests using Stochastic Weather Generator*; GeoConvention: Calgary, Canada, 2018.

36. Hosseini, S.M.; Ataie-Ashtiani, B.; Simmons, C.T. Spring hydrograph simulation of karstic aquifers: Impacts of variable recharge area, intermediate storage and memory effects. *J. Hydrol.* **2017**, *552*, 225–240. [CrossRef]

37. Dogančić, D.; Afrasiabian, A.; Kranjčić, N.; Đurin, B. Using Stable Isotope Analysis (δD and δ18O) and Tracing Tests to Characterize the Regional Hydrogeological Characteristics of Kazeroon County, Iran. *Water* **2020**, *12*, 2487. [CrossRef]

38. Ford, D.; Williams, P.D. *Karst Hydrogeology and Geomorphology*; John Wiley & Sons: West Sussex, UK, 2007.

39. Perne, M.; Covington, M.; Gabrovšek, F. Evolution of karst conduit networks in transition from pressurized flow to free-surface flow. *Hydrol. Earth Syst. Sci.* **2014**, *18*, 4617–4633. [CrossRef]

40. Zoghbi, C.; Basha, H. Simplified physically based models for free-surface flow in karst systems. *J. Hydrol.* **2019**, *578*, 124040. [CrossRef]

41. Mangin, A. Pour une meilleure connaissance des systèmes hydrologiques à partir des analyses corrélatoire et spectrale. *J. Hydrol.* **1984**, *67*, 25–43. [CrossRef]

42. Raju, K.S.; Kumar, D.N. Review of approaches for selection and ensembling of GCMs. *J. Water Clim. Chang.* **2020**, *11*, 577–599. [CrossRef]

43. Wang, H.M.; Chen, J.; Xu, C.Y.; Zhang, J.; Chen, H. A framework to quantify the uncertainty contribution of GCMs over multiple sources in hydrological impacts of climate change. *Earths Future* **2020**, *8*, e2020EF001602. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.