Hydrodynamic Characterization of Mugla Karst Aquifer Using Correlation and Spectral Analyses on the Rainfall and Springs Water-Level Time Series

Abstract

Karst aquifers have been an important research topic for hydrologists for years. Due to their high storage capacity, karst aquifers are an important source of water for the environment. On the other hand, it is safety-critical because of its role in floods. Mugla Karst Aquifer (SW, Turkey) is the only major water-bearing formation in the close environs of Mugla city. Flooding in the wet season occurs every year in the recharge plains. The aquifer discharges by the seaside springs in the Akyaka district which is the main touristic point of interest in the area. Non-porous irregular internal structures make the karsts more difficult to study. Therefore, many different methodologies have been developed over the years. In this study, unit hydrograph analysis, correlation and spectral analyses were applied on the rainfall and spring water-level time series data. Although advanced karst formations can be seen on the surface like the sinkholes, it has been revealed that the interior structure is not highly karstified. 100–130 days of regulation time was found. This shows that the Mugla Karst has quite inertial behavior. Yet, the storage of the aquifer system is quite high, and the late infiltration effect caused by alluvium plains was detected. This characterization of the hydrodynamic properties of the Mugla karst system represents an important step to consider the rational exploitation of its water resources in the near future.

1. Introduction

There are numerous methods to study, define and classify the karst aquifers. Statistical approaches are often in use. Some of them are well established and have been applied by the hydrogeologists for years. Time series analyses for forecasting and control were proposed by Box and Jenkins [1]. These statistical approaches were imported to karst environment researches by Mangin [2,3]. Then, they got widespread applications all over the world [4,5,6,7,8,9,10,11,12,13].

In the time series approach, karst aquifer system has an input and an output that correspond to the precipitation and discharge, respectively. The internal structure of the karst system acts as a black-box model. The black box modifies the input signal and creates an output signal. In the karst environment, this modification is mainly caused by the inner karstification level of the aquifer and by initial conditions of water level in aquifer. A high level of karstification creates a signal with less modification. Thus, karst development can be evaluated by analyzing the statistical relation between the precipitation and the discharge. Correlation and spectral analyses have been used to investigate input-output signal modification [14,15,16,17,18,19,20,21,22]. With the help of these statistical analyses, karst aquifers can be classified relatively with the other studied karst aquifers in the literature.

One quarter of Turkey’s land is covered by the carbonate rocks [23]. Anatolia hosts broadly distributed karst structures. Besides carbonate rock karsts, gypsum karst can be founded in a limited part of Turkey. Ekmekci [24] stated two major karst types in Turkey: evolutionary karst, which is predominant in Turkey, stands for continuous karstification with different maturation levels and rejuvenated karst, which implies reactivating a formerly developed karst.

Turkey has a long coastline cover. For this reason, karstification near the coasts was immensely affected by sea-level changes. Springs as discharge points of coastal karst aquifers are often located by the sea or in the sea. Bayari et al. [25] discovered 150 submarine groundwater discharges at the western coast of Antalya, SW Turkey. This situation threatens freshwater resources with salinization. The awareness for the protection and the effective management of freshwater has great importance in a country such as Turkey, which has a potential water scarcity in the future [26].

In this paper, Mugla karst aquifer is investigated in terms of karstification development level. Karstification level of Mugla karst gives crucial information about the aquifer in the area since it is the only groundwater resource in the region. In this context, time series analyses were applied on the Mugla karst to evaluate the hydrodynamic properties of the karst systems. This assessment will give a perspective on the karstification development level and contribute to a better understanding of the Mugla karst system. In line with this purpose, Mugla precipitation and springs water-level time series were used. Most of the published works on this matter deal with precipitation and spring discharge data. However, in some works, water-level data were used instead of discharge rates, due to the lack of the latter data [27]. In this study, no continuous discharge rate data of the springs were available. Thus, water-level was used. The correlation and spectral analyses were thus performed on the precipitation and water-level time series. In addition to time series analyses, recession coefficients for unit hydrographs were calculated to define the state of internal parts of the karst; matrix, fractures and conduits.

2. Materials and Methods

2.1. Study Area

The region is located within the Gökova Graben structure. The Gökova Graben region has been studied extensively in terms of geological evolution [28,29,30,31]. During the Miocene (23 million to 5.3 million years ago), the Gökova Graben region began to form and has been active since then [28]. Most studies focus on tectonic and active tectonics, sea-level changes, coastal types and physical geography of the region. Kurt et al. [29] discovered the buried Datça Fault located at the south and its antithetic faults located at the north of the Gulf of Gökova and discussed the relation between the basement unit (Lycian Nappes) and the Late Miocene-Plio-Quaternary basin fill. Gürer and Yılmaz [30] studied the geology of the close environ of the study area by indicating the local geological evolution. They proposed four episodes for the geological development. First, east (E)-west (W)-trending basin formed during Oligocene. Second, north (N)-south (S)-trending basin formed during Early Miocene. Third, a group of basins formed along north-northwest (NNW)-south-southeast (SSE)-trending faults during Late Miocene. Fourth, Gökova graben developed along E-W-trending normal faults. Gürer et al. [31] concluded that there are three basins formed by N-S compressional system and there are five basins formed by N-S extensional systems. Uluğ et al. [32] discussed five deltaic sequences at the northeastern Golf of Gökova and a relatively younger active faulting called the Gökova Transfer Fault.

The name of the karst spring group in the Gökova Bay is the Azmak springs. Akyaka, the neighborhood where the Azmak springs are located, is a part of the Gökova Bay. The Azmak spring group also creates a stream called Azmak. The 2 km long Azmak Stream discharges to the Gökova Bay, Aegean Sea (Figure 1). The Gökova Bay is designated as a special environmental protection area by the Republic of Turkey Ministry of Environment and Forests in 1988.

Although the altitude is quite variable in the province of Mugla, the elevation of the plains in the study area is between 600–650 m above the sea level (ASL) and the peaks vary between 800–900 m ASL. The Azmak Stream and the springs have an elevation between 0 and 10 m ASL. According to Turkish State Meteorological Service, there is a Mediterranean climate in areas between 0 and 800 m of elevation ASL, and Mediterranean mountain climate at higher elevations. The average annual rainfall is slightly over 1100 mm. Most of the precipitation falls in winter-spring and drought occurs in summer.

The Gulf of Gökova, where the Azmak Stream discharges, was formed as a result of a horst-graben system. The karst springs forming the Azmak Stream developed at the edge of the graben, the south of the horst. Akyaka, where the springs are in, is located on the hillside of a sharp elevation caused by the horst structure. The Azmak Stream is predominantly fed from karst aquifer, which has a 497 km² surface drainage basin (Figure 2). The Mugla city center and Mugla Plain (also called Karabaglar Plateau) region are not in the surface drainage basin of Azmak Stream. Karabaglar Plateau is a closed basin (Figure 2). In the literature, several studies investigated the hydrogeology of the area and the Azmak springs. Kurttaş [33] calculated the hydrological budget, revealing that there is additional input from the neighboring basin from the north (Mugla city center and Mugla Plain). In other words, the Mugla basin is in karstic connection with the Azmak springs. Kurttaş et al. [34] and Bayarı and Kurttaş [35] determined the undersea discharge points of the karstic springs on the coast of Gökova by using satellite remote sensing techniques and hydrological budget. Bayarı et al. [36] studied the recovery of freshwater discharges and the detection of cold-water points along the Mediterranean coast including the Gulf of Gökova by different techniques. Ekmekçi et al. [37] worked on the determination of seawater mixture in Gökova karstic springs by hydrochemical and stable isotope methods. They proposed that syphon and venturi processes are effective for saltwater intrusion for the Azmak springs. Acikel et al. [38] and Açıkel [39] formed the conceptual model of Gökova karstic springs using hydrogeology, hydrochemical and environmental isotopes. Additionally, in this conceptual model, the Mugla basin was included in the drainage basin of the springs. Acikel and Ekmekci [40] used multivariate statistical methods to evaluate groundwater quality of Azmak Springs. They suggested that the springs waters are divided into 3–4 groups created by the seawater intrusion effect and freshwater contribution. Yıldıztekin and Tuna [41] investigated the irrigation water quality of the well waters in the Mugla Basin and they found that the groundwater pollution was not a serious issue in the area. Sağir and Kurtuluş [42] demonstrated the performance of empirical Bayesian kriging and geo-adaptive neuro-fuzzy inference system on the Mugla Basin alluvium by revealing the covered/alluvium filled dolines in the area. Kurtuluş et al. [43] also studied on the groundwater of the Mugla Basin alluvium to assess the metal(loid) contamination. Yet, there is not any study evaluating the development level of the karst system.

The geological map of the area was digitalized after the MTA (General Directorate of Mineral Research and Exploration) geological map of Turkey (Figure 3) [44]. The most important geological units in which the karstification has developed are the Cretaceous-Late Triassic dolomitic limestone (JKmu) and the Early Jurassic-Late Triassic cherty dolomitic limestone (JKu). The JKmu outcrops at north and northwest of the area around Mugla, Akkaya and Yerkesik-Dogan plains while the JKu outcrops at the center and the east around Ula and Kizilagac plains. In Akyaka, mostly colluvium is present, and its contact with Middle Miocene aged post-tectonic conglomerate (Tk) can be seen. Additionally, JKu outcrops just over the springs in some places. Triassic-Permian schist (Pg) and Mesozoic peridotite (Kmo) outcrop at the center which has hydrogeological importance. Kmo is a part of Lycian nappes thrust that underlies JKu and overlies JKmu. The thickness of the Quaternary alluvium (Qal) in the plains is up to 100 m.

The most important factor of the region’s hydrogeology is the karstification of the carbonate rocks in the study area. The region shows most of the typical karst structures (Figure 4). Firstly, karst plains are important recharge areas in the region. The Mugla Plain is a closed basin, there is no stream coming out of the basin and all the water flows out of the basin as groundwater. Large and small karst structures are observed in all limestone-dolomite containing units in the region. The existence of sinkholes in all these plains is also known. The study region is quite active in tectonics, which creates faults, joints and fractures. These structural elements provided a favorable environment for karstification in the region where also rainfall is abundant. Since the horst-graben system in the region increases the hydraulic load difference between the recharge and discharge zones, it rules the groundwater flow and contributes to karstification. It can be suggested that the water drained in the Mugla Plain flows into the underground. Then, it is transferred through the southwest and Akkaya and Gulagzi Plains as groundwater flow. However, there is a low probability of groundwater flow from Mugla Plain to Ula Plain. The geological units can be cited as the reason for this. It is predicted that the Lycian Nappe thrust between the two plains prevents the groundwater flow. The Lycian Nappe unit, the Cretaceous peridotite and the serpentinized ophiolite are not suitable for groundwater flow. The Early Triassic metaconglomerate, sandstone, and siltstone unit (TRkv) in the same area is also impermeable. It is evaluated that the Lycian Nappe thrust continued towards the southwest of the region and was overlain by Middle Miocene post-tectonic conglomerate (Tk).

2.2. Data

Onset HOBO U20L-04 automatic data loggers were used for the water level measurements in the springs. This device can measure the water level between 0 and 4 m with ±0.4 cm accuracy and 0.14 cm resolution. The device can also measure the water temperature with an accuracy of ±0.44 °C and with a resolution of 0.1 °C. There are over 40 springs forming the Azmak Stream. However, the water level loggers were placed in L2 and L8 springs. L2 and L8 were chosen according to their temperature and electrical conductivity values measured over two years. These two springs showed different variations with time. Additionally, the rate of the discharges was enough during all the year when lots of the other springs dried out. The device placed at the L2 spring collected data between 02.09.2016–15.04.2018, and the device placed at the L8 spring collected data between 05.04.2017–05.02.2018. All the dates are given in dd.mm.yyyy format. All the measurements were done in 1-h intervals. The devices did not run out of battery during the measurement periods. However, since the memory of the devices got full after a certain number of measurements, the devices were removed from the springs and connected to the computer to transfer the data collected. Then, the loggers were placed back to the springs within no more than 20 min.

There are several weather stations near/in the study area. The rainfall data measured at the Mugla weather station represented the areal precipitation quite well. Açıkel [39] demonstrated that the regression is quite high (R² = 0.99) between the data of the Mugla station and the stations nearby. The Mugla station rainfall data (P) is used in the study (Turkish State Meteorological Service). The hourly precipitation data is from 00:00 01.01.2013 to 22:00 06.12.2017 (Figure 5).

Hourly water level (from the water surface to the logger) data were measured at the springs. “0 level” corresponds to the altitude ASL of the spring location. The altitudes of the L2 and L8 are 12 and 9 m ASL, respectively. The L2 water level data is from 10:00 02.09.2016 to 07:00 15.04.2018 (Figure 6a). The L8 water level data is from 16:00 05.04.2017 to 18:00 05.02.2018 (Figure 6b). Intersections of the data were used since all data are not exactly overlapped. Statistical analyses between L2 precipitation (L2_P) and L2 water level (L2_WL) were performed for the data between 10:00 02.09.2016 and 00:00 06.12.2017 with the data length of 10110 (Figure 7a). Statistical analyses between L8 precipitation (L8_P) and L8 water level (L8_WL) were performed for the data between 16:00 05.04.2017 and 00:00 06.12.2017 with the data length of 5509 (Figure 7b). General descriptive statistics and the box plot for L2_WL and L8_WL data are given in

Table 1. Descriptive statistics for L2 and L8 water level data.

(cm)	L2_WL	L8_WL	(cm)	L2_WL	L8_WL
Mean	26.91	24.86	Kurtosis	0.13	−0.26
Standard Error	0.05	0.08	Skewness	0.08	0.20
Median	26.55	24.76	Range	28.87	28.28
Mode	26.40	22.97	Minimum	11.11	11.67
Standard Deviation	4.71	5.81	Maximum	39.98	39.95
Sample Variance	22.18	33.75	Count	10110	5509

and Figure 8.

Distribution fitting and cumulative distribution fitting comparison for the L2_WL and L8_WL data were investigated in order to find out the statistical distribution characteristics and the similarities of the data. Distribution fittings were done by using maximum likelihood estimation with a 5% significance level. The L2_WL data were best fitted to Gamma-2 type distribution with the p-value of 3.08 × 10⁻¹⁴ (Figure 9a). The L8_WL data were best fitted to logistic type distribution with the p-value of 1.80 × 10⁻¹¹ (Figure 9b). The cumulative distributions of the L2_WL and the L8_WL data were compared by the two-sample Kolmogorov-Smirnov test with the significance level of 5% (Figure 10). The p-value was found lower than 0.0001, which indicates that the distributions of the data are not identical. Thus, using these two datasets to make different interpretations is significative.

2.3. Method

2.3.1. Recession Coefficient

All the calculations in this paper were carried out on Microsoft Excel with XLSTAT add-in.

Spring water level data were used for recession analyses to achieve a better understanding of the hydrodynamic properties of the karst aquifer. Since the data is hourly, monthly moving average values were used to demonstrate a clearer view of the data. Then, recession coefficient calculations performed using the Maillet [46] equation given in Equation (1):

$${\mathrm{Q}}_{\mathrm{t}}={\mathrm{Q}}_{R_0}{\times \mathrm{e}}^{-\mathsf{\alpha }\mathrm{t}}\to \mathsf{\alpha }=-\frac{{\mathrm{logQ}}_{\mathrm{t}}{-\mathrm{logQ}}_0}{\mathrm{t}\times \mathrm{loge}},$$

(1)

• α: Recession coefficient (1/day, multiplied by 24 for the conversion from hour to day)
• Q_t: Flow rate at the time t (Water-level in cm was used.)
• ${\mathrm{Q}}_{R_0}$: The last flow rate before the recession (Water-level in cm was used.)
• t: Recession duration (hours)

2.3.2. Simple/Autocorrelation and Simple Spectral Density Function

The simple correlation function is calculated by the following formula (Equation (2)) after Box et al. [47]:

$${\mathrm{r}}_{\mathrm{x}}=\frac{{\mathrm{C}}_{\mathrm{k}}}{{\mathrm{C}}_0}\to {\mathrm{C}}_{\mathrm{k}}={\mathrm{n}}^{-1}\sum _{\mathrm{i}=1}^{\mathrm{n}-\mathrm{k}}({\mathrm{x}}_{\mathrm{i}}-\overline{\mathrm{x}})({\mathrm{x}}_{1+\mathrm{k}}-\overline{\mathrm{x}}),$$

(2)

• r_k: Autocorrelation function
• k: Time lag, from 0 to m which is the cutting point
• C: Correlogram
• n: Length of the time series
• i: Time
• x: Parameter, water level
• $\overline{\mathrm{x}}$: Mean value of the time series

Spectral density function S(f) (Equation (3)), corresponds to the change from a time mode (time series space) to a frequency mode by applying the Fourier’s transformation on the variables [47]. The cut-off frequency is the time lag when the S(f) of water quantity data gets the values lower than 1 and gets negligible. This means that the input (precipitation) events before the calculated time lag are filtered.

$$\mathrm{S}\left(\mathrm{f}\right)=2\left[1+2{\displaystyle \sum }_{\mathrm{k}=1}^{\mathrm{m}}{\mathrm{D}}_{\mathrm{k}}{\mathrm{r}}_{\mathrm{k}}\cos \left(2\mathsf{\pi }\mathrm{Fk}\right)\right],$$

(3)

• S(f): Spectral density function
• F: is equal to j/2m
• j: is from 1 to m

D_k is a window being used to get an unbiased estimation of S(f). Mangin [2,3] suggested Tukey’s window as the best performing (Equation (4)):

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

(4)

The determining of the lag k and the truncation m is crucial for the correlation and the spectral analyses according to Mangin [2,3], who also suggested that it is much better to set an m value not higher than n/3, while n is the length of the time series. The value form is never selected to be higher than n/3 throughout this study.

2.3.3. Cross-Correlation and Cross-Spectral Density Functions

The cross-correlation function is calculated by the following formula after Larocque et al. [6] (Equation (5)):

$${\mathrm{r}}_{\mathrm{x},\mathrm{y}}\left(\mathrm{k}\right)=\frac{{\mathrm{C}}_{\mathrm{x},\mathrm{y}}\left(\mathrm{k}\right)}{{\mathsf{\sigma }}_{\mathrm{x}}{\mathsf{\sigma }}_{\mathrm{Y}}},{\mathrm{C}}_{\mathrm{x},\mathrm{y}}\left(\mathrm{k}\right)=\frac{1}{\mathrm{n}}\sum _{\mathrm{i}=1}^{\mathrm{n}-\mathrm{k}}({\mathrm{x}}_{\mathrm{i}}-\overline{\mathrm{x}})({\mathrm{y}}_{\mathrm{i}+\mathrm{k}}-\overline{\mathrm{y}}),$$

(5)

• r_x,y(k): Cross-correlation function
• C_x,y(k): Cross-correlogram
• σ_x: Standard deviation of the x time series
• σ_y: Standard deviation of the y time series
• $\overline{\mathrm{x}}$: Mean value of the x time series
• $\overline{\mathrm{y}}$: Mean value of the y time series

The higher degree of symmetry for the cross-correlogram graph indicates that there is another parameter affecting the variables. If the input signal is highly changed in the system, the maximum r_x,y(k) gets a lower value. In the karst environment, this means that the karst system is more complex and not highly developed. The lag to the maximum r_x,y(k) is called the delay which describes the transfer speed of the input change.

The cross-spectral density function, S_x,y(f) (Equations (6)–(8)), is the Fourier’s transformation of the cross-correlation function [6]:

$$\left|{\mathrm{S}}_{\mathrm{x},\mathrm{y}}\left(\mathrm{f}\right)\right|=\sqrt{{\mathrm{h}}_{\mathrm{x},\mathrm{y}}^2\left(\mathrm{f}\right){+\mathsf{\lambda }}_{\mathrm{x},\mathrm{y}}^2\left(\mathrm{f}\right)},$$

(6)

where:

$${\mathsf{\lambda }}_{\mathrm{x},\mathrm{y}}\left(\mathrm{f}\right)=2\left[{\mathrm{r}}_{\mathrm{x},\mathrm{y}}\left(0\right)+{\displaystyle \sum }_1^{\mathrm{m}}\left({\mathrm{r}}_{\mathrm{x},\mathrm{y}}\left(\mathrm{k}\right){+\mathrm{r}}_{\mathrm{y},\mathrm{x}}\left(\mathrm{k}\right)\right)\mathrm{D}\left(\mathrm{k}\right)\sin \left(2\mathsf{\pi }\mathrm{fk}\right)\right],$$

(7)

$${\mathrm{h}}_{\mathrm{x},\mathrm{y}}\left(\mathrm{f}\right)=2\left[{\mathrm{r}}_{\mathrm{x},\mathrm{y}}\left(0\right)+{\displaystyle \sum }_1^{\mathrm{m}}\left({\mathrm{r}}_{\mathrm{x},\mathrm{y}}\left(\mathrm{k}\right){+\mathrm{r}}_{\mathrm{y},\mathrm{x}}\left(\mathrm{k}\right)\right)\mathrm{D}\left(\mathrm{k}\right)\cos \left(2\mathsf{\pi }\mathrm{fk}\right)\right],$$

(8)

The coherence function CO_x,y(f) indicates the linearity of the input-output system (Equation (9)). The system is linear when the CO_x,y(f) near or equal to 1. High linearity is a signal of that any change in input has an analogical effect on the output. High linearity for karst systems’ input-output means well-developed karstification. Non-linearity shows there are other factors might be considered [6].

$${\mathrm{CO}}_{\mathrm{x},\mathrm{y}}\left(\mathrm{f}\right)=\frac{{\mathrm{S}}_{\mathrm{x},\mathrm{y}}\left(\mathrm{f}\right)}{\sqrt{{\mathrm{S}}_{\mathrm{x}}\left(\mathrm{f}\right){\mathrm{S}}_{\mathrm{y}}\left(\mathrm{f}\right)}},$$

(9)

If the gain function g_x,y(f) is higher than 1, there is an amplification of the output signal (Equation (10)). If the gain function g_x,y(f) gets a value lower than 1, there is an attenuation of the output signal [6]. The time lag for the amplification indicates the delay of the water release from the karst storage. While the discharge is being characterized by the precipitation, releasing the reservoir water beside the input related discharge creates an amplification after a certain delay (time lag of the gain function). The presence of amplification may also indicate the high storage capacity of the karst system.

$${\mathrm{g}}_{\mathrm{x},\mathrm{y}}\left(\mathrm{f}\right)=\frac{{\mathrm{S}}_{\mathrm{x},\mathrm{y}}\left(\mathrm{f}\right)}{{\mathrm{S}}_{\mathrm{x}}\left(\mathrm{f}\right)},$$

(10)

2.3.4. Karst Aquifer System Classification

Mangin [2,3] made a classification of karst aquifer systems in relation to the results of correlation and spectral analyses. This classification consists of four karsts that have been studied in detail before. This classification is based on four parameters: memory effect, spectral range, regulation time and unit hydrograph pattern (

Table 2. Karst aquifer classification proposed by Mangin [3].

Karst Types	Memory Effect (r = 0.1–0.2)	Spectral Range (Truncation Frequency)	Regulation Time	Shape of Unit Hydrograph
Aliou (Well Drained)	Poor (5 d)	Very Wide (0.30)	10–15 d
Baget	Small (10–15 d)	Wide (0.20)	20–30 d
Fontestorbes	Large (50–60 d)	Narrow (0.10)	50 d
Torcal (Poorly Drained)	Extensive (70 d)	Very Narrow (0.05)	70 d

).

The memory effect is the time lag for the correlogram reaches the value of 0.2. It points out the fall of the speed of the correlation. Rapid decrease exhibits the quasi-randomness of the response. The slow decrease means that the response is well structured. The slow decrease of the correlation implies a higher memory effect. The higher memory effect reveals that the karst’s idle nature. Idle karsts are often related to large storage. Hence, a correlogram decreasing slowly defines the poorly-karstified (poorly drained) system.

The regulation time corresponds to the duration of the input influence on the output. Regulation time is calculated by dividing the maximum spectral density value by 2, referring to the results of the spectral density function.

The spectral range corresponds to the frequency value when the spectral density is equal to zero. Broad spectral range points well-drained karst systems.

By using the mentioned approaches and methods above, a general view of the material, the method and expected results are given in

Table 3. Approaches, data, methods and the expected outputs.

Hydrograph Analysis		Autocorrelation				Cross-Correlation				Simple Spectral		Cross-Spectral
WL 1		P 2		WL		WL-WL		P-WL		WL		P-WL
Method	Output	Method	Output	Method	Output	Method	Output	Method	Output	Method	Output	Method	Output
Recession Coeff.	Differentiate karst zones, KL 3	Memory Effect	Dependence of data	Memory effect	KL	Symmetry check	Common input	Symmetry check	Detection of input-output relation	Spectral range	KL	Amplitude function, Cut-off frequency	KL
								Degree of the correlation	KL	Regulation Time	KL	Coherence function	Linearity of the system, KL
								Delay	KL	Cut-off frequency	KL	Gain function	Storage capacity, Late infiltration

¹ WL: Water level data. ² P: Precipitation data. ³ KL: Karstification level.

.

3. Results

3.1. Recession Coefficients

Two different recession coefficients for the L2_WL data were calculated in a selected period so that the water level is not disturbed by any rainfall (Figure 11a). The recession coefficients were calculated as 0.01941 day⁻¹ and 0.01174 day⁻¹ (Figure 11b). Two different coefficients indicate that the recession occurs in two phases. The first phase represents the respectively well-karstified zone (conduits/fractures) with higher transmissivity and hydraulic conductivity, and lower storage capacity. The second phase represents the respectively poorly karstified zone (matrix/fractures) with lower transmissivity and hydraulic conductivity, and higher storage capacity.

Three different recession coefficients for the L8_WL data were calculated in a selected period shown in Figure 12a. The recession coefficients were calculated as 0.01246 day⁻¹, 0.00956 day⁻¹ and 0.00854 day⁻¹ (Figure 12b). Three different coefficients indicate that the recession occurs in three phases. The first phase represents the respectively well-karstified zone (conduits) with higher transmissivity and hydraulic conductivity and lower storage capacity. The second phase represents the medium karstified zone (fractures) with medium transmissivity, hydraulic conductivity and storage capacity. The third phase represents the porous medium (matrix) with lower transmissivity and hydraulic conductivity and higher storage capacity. For this spring, the presence of the third phase means that the epikarst and the alluvium plains have an important effect on the water level of the spring discharge.

3.2. Simple/Autocorrelation

L2_P and L8_P autocorrelogram are given in Figure 13a,b. Autocorrelation r(k) gets a value of 0.2 at the time lag of 3 h for the L2_P data. r(k) gets 0.2 at 2 h for the L8_P data. The time lag for r(k) gets 0.2 represents the memory effect of the data. Too fast memory effect is an indication that the data is an independent variable showing quasi-randomness.

Autocorrelogram for the L2_WL and the L8_WL are given in Figure 14a,b. For the L2_WL, r(k) gets 0.2 at 1751 h, which makes 73 days as the memory effect. For the L8_WL, r(k) gets 0.2 at 1212 h which makes 51 days as the memory effect. Shorter the memory effect more developed the karstification. L2_WL memory effect is equal to 73 days, which is extensive and indicates poorly drained and karstified Torcal type karst aquifer (Table 2). The L8_WL memory effect of 51 days is large, which indicates better drained and karstified aquifer (Fontestorbes type) than the L2_WL. Moussu [48] listed the karst systems of Fontaine de Chartreux, Fontaine de Vaucluse and Durzon having the memory effects as 71, 75 and over 100 days, respectively. Lorette et al. [49] indicated 77 days of memory effect for Toulon spring.

3.3. Simple Spectral Density Function

The simple spectral density function for the L2_WL and the L8_WL are given in Figure 15a,b. S(f₀) values are 5233 and 6508 for the L2_WL and L8_WL respectively. For L2_WL data spectral range is lower than 0.1 and the regulation time is 109 days. For L8_WL data spectral range is again lower than 0.1 and the regulation time is 136 days. The cut-off frequencies are 3.1 and 3.5 days for the L2_WL and L8_WL respectively. This means that the rainfall takes place no longer than 3–3.5 days is being filtered by the karst aquifer system, and no effects of those rainfall events on the discharge can be detected. For both data, the spectral ranges are very narrow and the regulation times are too long, which are the indications of the Torcal type poorly drained and karstified aquifers. The classification points a much more poorly karstified aquifer than the Torcal since the parameters’ values are not very close to the Torcal type.

3.4. Cross-Correlation Function

The cross-correlogram between L2 and L8 water level data is given in Figure 16. The highly symmetrical plot of the cross-correlation for the L2 vs. L8 water level data indicates that there is a common input and the reaction against that input of these two data is similar.

The cross-correlograms of precipitation versus L2 and L8 water level data are given in Figure 17a,b respectively. Cross-correlation function for precipitation versus L2 and L8 discharge demonstrated a dissymmetrical plot, which is the mark of influence of the input (precipitation) on the output (discharge). The maximum cross-correlation values for L2 and L8 were calculated as 0.1564 and 0.1455, respectively. The low degree of the maximum r_x,y(k) values for the springs showed that the input was highly changed throughout the output. Slightly lower values for L8 might be the indication of a more complex recharge system.

“The delay” described as the transfer speed and corresponds to the time lag to the maximum r_x,y(k). The shorter delay indicates the faster transfer of the input influence on the output. The delays were 24 and 76 h for L8 and L2 respectively. To compare, Larocque et al. [6] calculated 0.5–0.8 for maximum r_x,y(k) values and 12–17 h of delay for a karst system, which had approximately 75 days of regulation time. Lo Russo et al. [50] calculated 1, 1, 1 and 38 days for different springs, which have 49, 57, 63 and 52 days of memory effect, respectively. Additionally, Fu et al. [51] calculated the 1-day delay for a discharge with a memory effect of 4 days. Thus, the signal transfer speed of the system is quite low which indicates a low degree of karstification.

3.5. Cross-Spectral Density Function

The amplitude functions (cross-spectral densities) for the L2_P-L2_WL and the L8_P-L8_WL are given in Figure 18a,b. Cut-off frequencies are 6.1 and 12.5 days for the L2_P-L2_WL and the L8_P-L8_WL. These values are quite high in comparison with Bange-L’Eau-Morte of 4.5 days [16] and Ain Zerga of 5 days [21]. This shows that the Mugla karst has quite inertial behavior.

The coherence functions for the L2_P-L2_WL and the L8_P-L8_WL are given in Figure 19a,b. The average coherence values are 0.21 for both L2 and L8. Low coherence means the non-linearity for the input-output system, and non-linearity might be a sign of poor karstification.

The gain functions for the L2_P-L2_WL and the L8_P-L8_WL are given in Figure 20a,b. For both data set, the gain functions are dominantly higher than 1 up to the 0.1 frequency. The amplifications occur almost at every frequency which shows that the karst system has high storage capacity because the system discharges the reserve water even at later frequencies. There are also the amplifications at mid frequencies indicating that there is very strong late infiltration.

4. Discussion and Conclusions

Correlation and spectral analyses were applied to the data of Mugla precipitation and the water level of two springs with the intent to achieve hydrodynamic properties of the karst aquifer. The karstification level/development of the karst and its internal parts were evaluated. The methods and the outputs are discussed as the same order as they were given in Table 3. All the expected outputs were achieved, even when the water level data were used instead of flow rate.

Hydrograph analysis were done using water level data of the two spring. The recession analysis in this paper was only used to discriminate the different number of phases for the springs if there are any. Since the used data are water levels of the springs, the calculated recession coefficients do not bear any physical sense compared to recession coefficients derived from flow rates analysis. The purpose was to show the different recession phases for the springs, which is important for understanding the system functioning. A complete recession analysis should split the flood flow and the base flow. However, in this paper, this is not the case, because there is no effective rain after the starting point of the recession. For that period, for L2, there are about five rain events with an average of 1–2 mm. For L8, it is more or less the same situation with a slightly higher amount of rain, of 2–3 mm. Taking into account the fact that the numerical values of the recession coefficients are not significant, the rainfall amount in the period of recession is already negligible to consider the infiltration effect. According to Mangin [52], the infiltration effect can be neglected under these kinds of conditions. In order to calculate the dynamic volume, storage and hydraulic conductivity through the recession curve, flow rate data are required. Thus, it was not possible to calculate these parameters using conventional methods with water level data. Horoi [27] also used the water level data instead of discharge rate data of the springs for the statistical analyses. Since most of the calculations are unitless, this may not cause any error. However, the relation between the water level and the flow rate is not proportional. This issue should be addressed by further investigations comparing these two kinds of data for the same discharge. Considering the recession analysis, karst internal parts were grouped as matrix (alluvium plains/dolines/poljes), fractures (poorly karstified zone) and conduits (well-karstified zone). For the L2 spring, only matrix and fracture parts of the karst were determined. For the L8 spring, all three main parts of the karst were detected. The calculated recession coefficients might indicate poorly developed conduit structure and a low degree of karstification.

The uncertainty caused by the use of the water level data instead of the flow rates is also valid for the correlation and the spectral analyses results. However, in these analyses, the calculations and their values bear physical sense. The calculated values are in accordance with those in the literature. Autocorrelation and simple spectral calculations using water level data provided memory effect, spectral range, regulation time and cut-off frequencies. Cross-correlation and cross-spectral calculations using the precipitation and the water level data provided delay, degree of correlation, linearity of the system, storage capacity and late infiltration effect. Water level data were used for all these calculations also. Yet, the results were consistent with each other and those in the literature and had hydrogeological meanings, which indicated a low degree of karstification for the studied karst system. According to the correlation and spectral analyses, the response (spring discharge) of the input (precipitation) is slow, non-linear and well modified, which are again the signs of poor karstification. Yet, the karst aquifer has high storage capacity and a strong late infiltration effect. The important late infiltration effect stands for the dominant role of the matrix zone.

In light of the findings on the Mugla karst system, it can be suggested that the karst development level of this system is low. It can be concluded that the epikarst and the alluvial plains are still the most important parts of this karst system. The Mugla karst system is a vital water resource for the region. It is therefore essential to pay more attention to these parts of the karst system while managing it in a sustainable way. Considering the revealed features of the karst system, it is not possible to create sudden floods in the discharge area, but attention should be paid to the flooding that may occur due to late drainage of water in the feeding area.