Characterization of the hydrodynamic functioning of the Degracias-Sicó Karst Aquifer, Portugal

A study on the complex processes of hydrogeological transformation of recharge into discharge within a karst system is described. Given the absence of boreholes or points of direct observation, the applied methodology is based on study of the natural response of the main outlet of the karst aquifer, considered here as a black-box. Time series analysis and recession curve analysis are applied to daily time series of rainfall as input, and spring discharge as the output of the aquifer. Daily variations of temperature and electrical conductivity of spring water are considered as well. The present study, an integrated research endeavor, is the first quantitative approach toward understanding the hydrogeological functioning of a large and not-well-known karst aquifer and improving existing knowledge about its internal physical characteristics. The case study is the Degracias-Sicó karst aquifer, in Portugal’s central-western region, which is relatively unknown from the hydrogeological perspective. The response of the karst aquifer’s main spring shows an important component of quick flow after heavy rainfall events and then a very slow recession. This bimodal behavior is proof of the significant spatial complexity and heterogeneity of the karst aquifer, marked by the presence of a conduit-dominated flow throughout a connected conduit network system and the dominance of base flow that confirms the aquifer’s large storage capacity. The results also demonstrate a significant component of delayed flow.


Introduction
Karst aquifers are typically characterized by a high heterogeneity due to the irregular network of pores, fissures, fractures and conduits with very different forms, sizes and connections, which make them unique and different from other aquifers (Bakalowicz 2008). Such structural uniqueness gives rise to a very specific and complex hydrodynamic functioning and to hydrogeological parameters that can vary enormously over space and time. The specificity of karst aquifers is also on the duality of recharge, storage and circulation processes (Smart and Hobbs 1986). Infiltration into the karst system can be concentrated via sinkholes (swallow holes) and boreholes and/or diffused through fractures and fissures in the outcropping rock or the soil. This anisotropy is responsible for different groundwater flow paths and velocities. The first occurs in the conduits and enlarged fractures networks, which allow a quick transit through the vadose zone to the saturated zone. The diffuse infiltration starts a slow percolation through the small fissures in the vadose zone. According to its thickness and specific conductibility characteristics, the infiltrating water can take a long time to reach the saturated zone.
Based on the level of development of the conduit network and its interconnectedness within the karst system, it is possible to define two types of karst aquifers (White 1988;Ford and Williams 2007;Taylor and Greene 2008). One is the conduit karst system characterized by a well-developed and wellorganized network of karst conduits, which allows for quick drainage after a rainfall event due to the high transmissivity and low storage capacity. The result is an intense discharge of the karst spring over a short period represented in the hydrograph by a peak discharge. These aquifers are known as quick flow or fast flow karst systems and correspond to the typical karst aquifer definition. In contrast, the diffuse karst system differs from the first type due to the dominance of base flow at the expense of fast flow. In this type of karst system, the conduit network is poorly interconnected or not well developed, so that water takes a long time to travel from the recharge area to the outlet. These diffuse karst aquifers retain water for a much more extended period and have a large capacity to store water, and to release it later during a dry period or dry season.
The heterogeneity, duality, anisotropy, nonlinearity and nonstationarity that characterize karst aquifers make their study a real challenge (White 2002). Hydrogeological investigation into karst aquifers is particularly hard due to the complex nature of karst systems and the limitations or difficulty of applying Darcy's law and classical hydrogeological investigation techniques such as well observations and measurements, pumping tests and tracer tests (White 1988;Ford and Williams 2007;Goldscheider and Drew 2007;Andreo et al. 2008;Kresic 2010;Fiorillo 2014;Mikszewski and Kresic 2015). These intrinsic complexities require the use of specific methods (Goldscheider and Drew 2007;Taylor and Greene 2008) to provide information on the different types of flow, groundwater circulation pathways, storage capacity and spatial and temporal variations. Thus, indirect methods of characterization have been applied to draw maximum information from the karst systems (black box model, used in this study, and distributed parameters approach). Relevant reviews and detailed discussions about the methods to investigate karst hydrodynamics are presented for instance in Bakalowicz 2005;Kovács et al. 2005;Ford and Williams 2007;Goldscheider and Drew 2007;Taylor and Greene 2008;Martínez-Santos and Andreu 2010;Rehrl and Birk 2010;Kresic and Stevanovic 2010;Ghasemizadeh et al. 2012;Geyer et al. 2013;Hartmann et al. 2014).
Investigation techniques for understanding karst hydrodynamics, used in preliminary studies performed to acquire basic information, are based mainly on the analysis of output signals at the system outlet that reflect the karst system's specific settings (Grasso and Jeannin 2002). The spring discharge, as the global system response, is the combination of the different processes associated with recharge, storage, and throughflow in a karst basin upstream from its outlet (Kresic 1997). Therefore, water discharging from a karst spring carries an imprint of the global hydrodynamic functioning of the karst aquifer and consequently, of its physical features.
Focusing on the analysis of a spring discharge, this study applies a combined methodological approach that is widely used in research on the hydrogeological functioning of karst hydrosystems, which including time series analysis on rainfall and discharge, recession curve analysis and physical parameters analysis on spring discharge. The present investigation aims to provide in-depth knowledge, in functional terms, of the hydrodynamic behavior of an almost unknown karst aquifer, in order to understand the transfer processes occurring within the karst system from recharge to discharge. Moreover, this study intends to infer information about the functionality of the underground karst aquifer's network drainage system, and its storage capacity as well as the internal degree of karstification.

Study area
The Degracias-Sicó karst aquifer is located in the centralwestern region of Portugal, 30 km from the Atlantic Ocean and a few kilometers south of Coimbra, covering about 100 km 2 in the northern part of the Sicó karst massif (Fig. 1). A set of limestone mountains and a central plateau shape the morphology of this area. The highest points are at the Serra do Rabaçal mountain (532 m above sea level (m asl) at the northern border, and at the Serra do Sicó mountain (555 m asl) at the southern border. The central sector of the massif is dominated by Degracias-Alvorge plateau with elevations ranging from 280 to 350 m asl. This roughly flat surface corresponds to a wide morpho-structural depression that favors the infiltration process of the endorheic basin.
Geologically, the bedrock of the northern area of the Sicó massif is composed of a thick well-stratified and generally continuous Jurassic limestone succession organized in a large monoclinal tectonic structure dipping towards the west (Cunha 1990). The post-Jurassic tectonic deformations were responsible for the progressive uplift of the Sicó massif, as well as for its consequent internal compartmentalization (intense block-faulting). The western border has been affected by important fault systems, mainly N-S to NNE-SSW and NE-SW to ENE-WSW directions (Fig. 1), which play an essential role in the groundwater circulation in the massif.
Compact limestones, oolitic limestones and dolomiticlimestones of Degracias and Senhora da Estrela formations (early Bajocian-Bathonian) are the lithostratigraphic units where the surface/subsurface karst features are more abundant and pronounced in the Sicó massif (Soares et al. 2007;Cunha and Dimuccio 2017). In contrast with the uncovered carbonate rock, most of the Degracias-Alvorge plateau (in central and eastern areas) is covered by remains of the siliciclastic cover (Lower Cretaceous), which discordantly overlaying Jurassic carbonates.
The climate in the region is the temperate Mediterranean. The mean annual rainfall recorded from 1979 to 2008 at the precipitation station inside the study area is 1,200 mm. This amount is corroborated by the study's meteorological station (installed in 2009 at a very close location) and demonstrates the orographic effect (relief barrier) of the Sicó karst massif. The maximum rainfall occurs in the November-March period (57.8% of the total annual), resulting from the latitudinal atmosphere instability and wet winds from the Atlantic Ocean, and the minimum rainfall is in July−August (3.6%). The wettest trimester (November, December, and January) amounts to 41.3% of the annual rainfall, while the driest trimester (June-August) does not exceed 7.5%. The real evapotranspiration computed using different methods- Turc (1954), Coutagne (1954), Thornthwaite and Mather (1957), Kessler (1965)-is higher than 50% of annual rainfall and is opposite to the temporal pattern of rainfall. During the spring-summer season, the rate of evapotranspiration is so high that the recharge can be considered null, whereas in the autumn/winter season, due to the very low rate of evapotranspiration, the monthly recharge may reach values higher than 82% of rainfall .
The spatial distribution pattern of infiltration water is substantially different in terms of quantity and velocity. The aquifer recharge is more concentrated and quicker in the wellkarstified areas and uncovered rock of Serra de Sicó and in the faulted western border of Degracias Alvorge plateau. In the areas covered by muddy-sand and conglomerate siliciclastics deposits, the soils exist up to a few meters thick, which is responsible for a diffuse and much slower infiltration. Recent studies carried out in the Degracias-Alvorge plateau (Paiva 2015) have proved that the recharge is totally autogenic and mostly diffuse, flowing through the fissures of the carbonate rocks and essentially through its karst landforms. Only locally, the recharge is concentrated, especially during significant rainfall events, when part of the water infiltrates immediately through highly karstified limestones. The endorheic features and the elevation of this area (300-550 m asl) and the surrounding areas explain the purely autogenic character of recharge.
Regarding the hydrodynamic functioning, the Sicó massif reveals an important structural control, particularly through its monoclinal geological structure, as well as through the fracturing system, which influences the regional groundwater flow towards the western border. Consequently, the largest springs (perennial and temporary) are located along the base of the western boundary of the massif (ranging from 40 to 75 m asl), establishing an alignment that follows the submeridian orientation of the fracturing system (Fig. 1).
The discharge occurs mainly through two perennial springs: Olhos d'Água do Anços (OAA spring), which is the main outlet of the hydrosystem, and Ourão Spring. According to the long-term records, both springs drain jointly from 75 to 80% of the water circulating in the massif in normal and wet years. During the dry years, that value frequently exceeds 90%. In the study period (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016), the mean annual volume drained by these springs was 77 Mm 3 (43.0 Mm 3 by OAA spring and 34.0 Mm 3 by Ourão Spring). The remaining springs are very small (temporary upstream to the main spring and perennial downstream) and drained together, accounting for only about 7 Mm 3 . These secondary outlets present a similar hydrodynamic behavior to the main spring, particularly after intense rainfall episodes.
The mean annual discharge of OAA spring is 1.4 m 3 /s, and the mean daily flow oscillates from less than 0.2 m 3 /s in summer to 5.8 m 3 /s in winter, which reflects the annual rainfall distribution (Fig. 3). The mean annual outflow of Ourão spring is 0.6 m 3 /s and the mean daily discharge is 0.2 and 2.7 m 3 /s in the summer and winter, respectively. Although it is important for the aquifer drainage, the Ourão spring was not considered in the detailed analysis, since its natural behavior has been altered by the construction of a water-retaining structure immediately downstream from the spring. Because of this, the natural responses of the karst system are unreliable and compromise the input-output analysis.

Data and methodology
A gauge station was installed at the main spring of the aquifer (OAA spring) to monitor parameters such as spring flow level (meters), temperature (°C) and electrical conductivity (μS/cm) of spring water. The data were collected by a Level-Troll 500 datalogger and a CTD Diver datalogger with an acquisition time interval of 20 min, from 1 October 2009 to 30 September 2016. The accuracy of the temperature sensor is ±0.1°C and that of the electrical conductivity sensor is ±1% of reading, with a minimum of 10 μS/cm. The second largest outlet of the Degracias-Sicó aquifer (Ourão spring) was not considered, as explained before, and the minor springs of the aquifer's western border are very small or intermittent, making it impossible to install a gauge station. The aim of this detailed information is to identify with precision the temporal variations of the spring water parameters. Rainfall has been recorded to 0.2 mm accuracy by a single rain gauge (Sra. da Estrela) situated at 350 m asl inside the study area, which can be considered representative of the rainfall events in the area, considering the region's air circulation pattern, topographic features and results from previous studies (e.g., Paiva 2015). The rainfall data have been validated by comparing the datasets to rainfall time series from rain gauges placed nearby and under similar geographical conditions; the rainfall data are available from the Portuguese Environment Agency's meteorological network. The location of the meteorological station is given in Fig. 1. The rainfall measured by this rain gauge is really demonstrative of the water available for the recharge of the karst system.
Piezometric level data and other useful hydrogeological information usually taken from wells, natural holes or groundwater conduit systems proved impossible to collect in this case, due to the physical features of the karst aquifer under study. The huge thickness of the karstic massif (more than 300 m between the Degracias-Alvorge plateau and springs level, in the massif's western border) makes the installation of piezometers impossible. By the same token, it is not possible to collect information about water circulation in the unsaturated zone or the fluctuations in the piezometric surface.
In the present study, the Degracias-Sicó karst aquifer is considered as a black box, by reason of the scarcity of information on the speleological characteristics of the Sicó karst massif, the very narrow geographic scope of previous scientific studies, and the inexistence of data series (rainfall and spring outflow). In order to study the karst aquifer hydrodynamics and to achieve some in-depth knowledge of its physical features, an input-output approach has been used. Methodologically, total rainfall was used without taking evapotranspiration into account. This was because of: (1) the scarcity of deep soil and dense vegetation in the recharge area (except for the eastern area, which has a significant soil layer); (2) the temporal coincidence between the highest amount of rainfall and the lowest temperature values (from November to March); (3) the flat nature of the recharge area, which tends to amplify the infiltration process; (4) the inexistence of relevant water bodies on the surface. Thus, the climatic and physical features of the recharge area are favorable to the infiltration process rather than cause evapotranspiration losses.
As the outlet of the karst aquifer, the spring discharge functions are an integrated representation of its overall behavior. The responses of the springs on the western shore of Sicó karst massif were taken into consideration, particularly in the case of the main spring. To arrive at a preliminary understanding of the functioning of the Degracias-Sicó karst aquifer, mean daily spring discharge was analyzed using three different methods: time series analysis (correlation and spectral analysis), spring flow recession curve analysis, and spring water temperature and electrical conductivity analysis. The use of an indirect methodology based on the spring discharge data is justified in the present investigation by a number of reasons: (1) there is no information available about the internal characteristics of the aquifer and data series; (2) it is easy to implement (data acquisition by rain and spring gauge stations installed at the start of the investigation); (3) it yields quantitative and reliable results, and (4) it has been widely tested. A short overview of the statistical methods used in this study is presented in the following subsections to make the results easier to understand and interpret.

Time series analysis
Time series analyses are signal processing methods and are commonly performed as preliminary studies to obtain a basic understanding of karst aquifers. Based on a systemic approach, this methodology was introduced by Mangin (1975Mangin ( , 1984 to research input-output relations in karst aquifers. Correlation (time domain) and spectral analyses (frequency scale) are one type of time series analysis for studying transfer functions between input (rainfall) and output (spring discharge) using statistical functions, including autocorrelation, spectral density, cross-correlation, cross-amplitude, gain, coherence and phase functions. For detailed theoretical explanations and mathematical expressions, the reader is referred to Jenkins and Watts (1968), Box et al. (1994), Mangin (1975Mangin ( , 1984, Padillla and Pulido-Bosch (1995), Larocque et al. (1998) and Massei et al. (2006).

Univariate analysis
The univariate analysis characterizes the individual structure of time series-autocorrelation in the time domain, and spectral density in the frequency domain. The autocorrelation function allows the quantification of the memory effect of a system. This function is represented by means of a correlogram and is computed by the decorrelation lag time up to the predetermined value of r (k) = 0.2 (Mangin 1975). The slope of the correlogram is determined by the system's response to a given event, and the quicker and sharper the decrease in the correlogram, the shorter the influence of the event on the karst system's response. The spectral density function is complementary to correlation analysis and corresponds to a change from a time mode to a frequency mode through a Fourier transformation of the autocorrelation function. The various peaks in the graph represent periodical phenomena in the output series.

Bivariate analysis
Bivariate analysis centers on the transformation of the input signal into an output signal. The cross-correlation function is examined in a time domain, whereas the cross-amplitude, phase, coherence and gain functions are analyzed in a frequency domain. The cross-correlation between precipitation and karst spring discharge gives an idea of the response to a single event (Padillla and Pulido-Bosch 1995;Larocque et al. 1998) and consequently provides an estimation of the degree of karstification of the aquifer (Panagopoulos and Lambrakis 2006).
The cross-amplitude function indicates the system's filtering features for the periodic components of the rainfall data (Padillla and Pulido-Bosch 1995). The karst systems that attenuate most of the short-term input signals (high frequencies) and amplify the long-term signals (low frequencies), are usually considered poorly developed karst systems, with high storage capacity and a nonefficient drainage network. The coherence function highlights the linearity of the karst system, which can be considered linear in that a change in the input produces a proportional change in the output (Larocque et al. 1998). The phase function can be used to estimate the mean delay, for different frequencies, between precipitation and spring discharge (Padillla and Pulido-Bosch 1995). The truncation point used for calculating all correlation and spectral functions was 100 days (see Rahnemaei et al. 2005).

Recession curve analysis
Hydrological time series analysis has proved to be especially useful in the study of karst environments, but it alone does not lead to a satisfactory understanding of the hydrodynamic behavior and characteristics of karst systems (Jeannin and Sauter 1998;Labat et al. 2000;Gárfias-Soliz et al. 2010).
Among the methods used to study the recession curve of karst springs, the classification method proposed by Mangin (1975) is the most achieved because it characterizes the infiltration conditions quantitatively throughout the unsaturated zone and the storage capacity of the saturated zone (Marsaud 1997;El-Hakim and Bakalowicz 2007). According Mangin (1975), a nonlinear flood recession, expressed by the function ψ (t) , represents the discharge from the nonsaturated zone; an exponential base flow recession, stated by the function φ (t) , expresses the discharge from the phreatic zone (Fig. 2). Thus, the entire recession curve can be presented according to Eq. (1): Q (t) is the discharge at time t, ψ (t) corresponds to the quickflow part of the recession and φ (t) to the base flow segment (Padilla et al. 1994). The base flow recession follows Maillet's expression (Ford and Williams 2007) as shown in Fig. 2 and expressed here as Eq. (2): where Q R0 is the discharge at the beginning of the slow recession and is extrapolated from t i (when the flood recession assumed essentially zero, i.e. ψ (t) = 0), and α is the base flow coefficient in day −1 . In order to characterize the infiltration, Mangin (1975) suggests using homographic function ψ (t) expressed by Eq.(3) (see Fig. 2): in which q 0 is the difference between the peak discharge (Q 0 ) and the discharge at the beginning of the slow recession (Q R0 ). Additionally, η = 1/t i given in day −1 , with [0, t i ] equal to the duration of the fast recession; η is related to the velocity of infiltration, thus the lower the value of η, the slower the infiltration. ε is the coefficient of heterogeneity, characterizing the concavity of the hydrograph and thus the type of groundwater circulation. The total initial volume (V 0 ) provided by the spring at t 0 (flood peak) is the sum of the dynamic volume (V dyn ) and the volume drained during the fast recession ðV * 1 ).
The dynamic volume V dyn under base flow conditions is calculated by integrating Eq. (4) (see Fig. 2).
The total volume drained by the spring during the fast recession is given in Eq. (5), also displayed in Fig. 2. Mangin (1975Mangin ( , 1984 also introduced two indices to classify the karst systems. k defines the extent and the regulation power of the saturated zone and implicitly indicates the mean residence time of water inside the karst system. i characterizes the infiltration conditions and represents the delay between recharge and spring outflow. The indices k and i are calculated by the Eqs. (6) and (7), respectively: Plotting i against k, Mangin (1984) classified karst aquifers into five groups. To determine the necessary parameters to implement Mangin's method, the Curve Fitting Toolbox from MATLAB software has been used. The recession curve analysis has been applied to daily spring recharge in seven seasonal recession periods, from 2010 to 2016.

Water temperature (T) and electrical conductivity (EC) analysis
The analysis of variations in water temperature and electrical conductivity at the karst springs permits characterization of the karst aquifer's behavior, as these parameters are nonconservative tracers. The temporal oscillation in these parameters permits one to differentiate immediately the drainage patterns, as conduit flow displays significant variations and diffuse flow exhibits low variations. The flow throughout the conduits system is considered 'new water', freshly infiltrated water, characterized by low electrical conductivity and temperature that differ from the thermal characteristics of water drained at that moment by the spring. Conversely, the 'old water' has a longer stay within the aquifer and often presents much higher temperature and electrical conductivity values due to the water-rock interaction.
In this study, water temperature and electrical conductivity were monitored from 1 October 2011 to 17 July 2014, at the OAA spring, with a resolution of 20 min using a Baro Diver, TD Diver and CTD Diver (sensitivity: 1 μS/cm and ± 0.1°C). The conductivity meter considered automatically the electrical conductivity values as corresponding to a standard temperature of 25°C. In order to assess the total reliability of the data series and to detect any possible anomaly of the divers, data of temperature and electrical conductivity were collected every 2 weeks, during the same period, in all the perennial and temporary small springs at the west border of the aquifer. The equipment was a Hanna multi-parameter portable meter Fig. 2 Idealized diagram representing the components of a karst spring recession curve according to Mangin's analysis method (after Ford and Williams 2007). Shaded area signifies base flow; unshaded area represents the quick flow. See text for the explanation of each formula and identification of variables (pH/T/TDS/EC), also calibrated to 25°C. Unfortunately, it was not possible to collect and analyze enough samples of spring water under different hydrometeorological conditions to determine important hydrochemical parameters, which would have brought valuable information to the analysis of the system's hydrodynamics.

Results and interpretation
Olhos d'Água do Anços spring hydrograph The daily hydrograph of OAA spring is characterized by pronounced peaks of discharge following heavy rainfall events, as shown in Fig. 3. The intensity of the spring response to a rainfall event is linked to the time the event occurs in the hydrologic cycle, that is, the cumulative rainfall so far and the hydraulic conditions of the aquifer. The rising is fast and attests to the enormous increase of the spring outflow, particularly in the case that the storm occurs when the underground conduit system is empty. In these situations, the increase of the spring flow begins on the same day of the recharge event (a short lag time), and the total outflow becomes several times higher. The hydrograph of OAA spring represents a karst-type response and indicates developed karstification inside the aquifer with a good connection of the conduit networks with the spring and the presence of concentrated recharge through active shafts. This assumption is reinforced by similar hydrodynamics observed in the small springs upstream and downstream of the OAA spring. The rising of spring flow is much quicker than the falling, which can be considered a characteristic of the aquifer dynamics as occurs independently of the hydrometeorological features of the year and the hydraulic conditions of the aquifer.
The recession curve presents a less steep limb than the rising limb for short periods and particularly over the long seasonal recession, which occurs during the late spring-summer period, extending into the following hydrological year. The smoother decrease of the spring flow from June and its relative stability during the dry season demonstrate the enormous storage capacity of the aquifer. The base flow which exclusively supplies the spring during the low-flow period corresponds to the water stored throughout the wet season and gradually released by the aquifer. The single rain events have no substantial influence on the spring discharge, especially if they occur from June to the end of September. The high temperature increases the evapotranspiration to the maximum, and infiltration becomes null.

Time series analysis
The autocorrelation function (ACF) of daily precipitation diminishes very quickly and reaches the confidence level (r k = 0.2) after 2 days. This means that no memory effect is detected in daily rainfall, which can be considered as a pure stochastic and unstructured phenomenon (pure aleatory function). The daily discharge series of OAA spring reaches r k = 0.2 after 79 days (Fig. 4). This shows the high memory effect of the system, according to the theoretical values (Mangin 1984) and the studies that use multiple-year time series (e.g., 83 days, Panagopoulos and Lambrakis 2006;76 days, Larocque et al. 1998;73 days, Jemcov and Petric 2009;66 days, Amraoui et al. 2003). The general slow decrease in the correlogram means that storage is significant, and water is released from the aquifer in a gradual manner. Furthermore, the correlogram Fig. 3 The hydrograph of daily mean discharge (m 3 /s) of Olhos d'Água do Anços spring with the total daily rainfall (mm) at Sra. da Estrela rain gauge during seven hydrological years (from 1 October 2009 to 30 September 2016) indicates a quicker drop at first (for less than 14 days), followed by a much slower decrease. The first section seems to correspond to the influence of karst conduits and enlarged fractures, which rapidly drain the aquifer. This inference seems to be corroborated by the correlogram behavior in dry years (e.g., 2011/2012), when the conduit system is empty and the water from a heavy rainfall event travels quickly through the unsaturated zone. The autocorrelation completely disappears (r k = 0.0) beyond the 100-day period, which could be attributed to the influence of aquifer reserves and very narrow conduits and discontinuities, which slowly empty the aquifer over a significantly long period.
The simple spectral density function (S f ) for the daily discharge of the OAA spring shows a striking peak at the frequency 0.0023 (432 days), which indicates the presence of an important cycle over more than 1 year. This can be partially explained by the annual recharge cycle of the aquifer (large peak rainfall daily data at S f = 0.00277, i.e. 360 days) and shows the significant pluriannual regulation of the reserves (Mangin 1984), which in turn confirms the huge storage capacity of the aquifer. The cut-off frequency is large (S f < 1 = 0.197, i.e. 5 days) which is not in accordance with the characteristics presented previously. A possible interpretation could be that a significant quick flow component, caused by the fast circulation in a connected conduit system, is superimposed on a dominant base flow component.
The cross-correlation function (CCF) for the spring daily data series displays the aquifer's bimodal hydrodynamic behavior expressed in a sharp peak, followed by a gradual decrease (Fig. 5a). This shape confirms the aquifer's duality, with a transmissive function being indicated by the sharp peak and a capacitive function represented by a slow decrease. The peak indicates a significant impulse response by the aquifer to the recharge and an incomplete filtration of the input signal. The shorter response time indicates a well-developed karst system, with a connected underground conduit network. After the peak, CCF decreases very slowly, which proves the aquifer's large storage capacity and its progressive emptying, as well as the dominance of the base flow component in the spring outflow. The time lag determined from the crosscorrelation function is 3.5 days, which is an indication of the aquifer's rapid response to a rainfall event and the presence of a significant quick flow component in the spring flow. The clear dissymmetry of cross-correlation values towards the positive k values confirms that rainfall unequivocally influences the spring's flow rates. The duality of the hydrodynamic functioning of the Degracias-Sicó karst aquifer is also made plain by the results of the cross-amplitude function (CAF; Fig. 5b). In the frequency domain, the CAF is characterized by a significant attenuation of the input signal as the maximum value occurs at low frequencies (0.0023, i.e. 432 days). This peak clearly indicates a cycle of more than 1 year, which emphasizes the aquifer's enormous storage capacity. Additionally, CAF has a large spectral band, as it only reaches values close to zero at frequencies higher than 0.382 (3 days) which in turn confirms the preceding statement and indicates the remarkable influence of the quick flow component on the spring's total outflow.
The coherence function (COF) expresses the linearity of the input-output relationship. The mean value of the daily rainfall-discharge coherence function (0.54) once again shows the complex hydrodynamic functioning of the aquifer. This value could be interpreted as the aquifer having some linearity in the rainfall-spring flow relationship and acting as a filter that attenuates the input signal as well. Significantly high values of COF for high frequencies (Fig. 5c) indicate that the aquifer has a linear behavior in response to the input. At the same time, the COF has considerably low values for low frequencies, reflecting a very well-regulated system (Panagopoulos and Lambrakis 2009). Finally, the phase function (PhF) shows a satisfactory alignment for low frequencies and a high distortion for high frequencies (Fig. 5d), which indicates the predominance of the base flow in the spring outflow. For frequencies below 0.07 (14 days), a better coherence can be observed. This is visible and acceptable for frequencies until 0.28 (3.5 days), which proves that the mean delay in the input-output signals is 3-4 days.
The spectral analysis of daily data series of rainfall and spring flow has clearly shown the presence of secondary peaks at frequencies that correspond to 17-, 39-, 57-and 90-day periods. Similar temporal values were identified on the correlograms of auto and cross-correlation. These temporal markers seem to indicate the existence of several welldefined responses that correspond to the arrival to the spring of flow with different travel times and characteristics. The explanation to this delayed flow could be: (1) the predominance of diffuse recharge; (2) the emptying of different sectors of the aquifer by the establishment of normal gradient from the matrix and discontinuities to the conduits that become unsaturated; (3) the contribution from the epikarst that covers the eastern and central parts of the recharge area; (4) the dimension of the aquifer which dictates the different travel times.

Recession curve analysis
The application of Mangin's method strengthened the previous results, which emphasizes the duality of the hydrological functioning of the aquifer. The results for the ψ (t) function indicate the prevalence of diffuse recharge and slow circulation in the unsaturated zone (η = 0.023; ε = 0.143) especially in years in which the rainfall amount is higher than the mean value (η = 0.010, ε = 0.033 in 2013; η = 0.012, ε = 0.032 in 2014; Table 1). The mean value of the i function (0.78) corroborates this statement, as i ranges between 0 and 1, and a result closer to 1 represents a karst system recharged by a slow or delayed infiltration. However, a significant quick flow component also is present in OAA spring discharge, expressed by the relevant values of ε, as the more concave the hydrograph, the more dominant the fast infiltration. This type of circulation is particularly observed in dry years. The concavity of the fast recession limb is ten times larger in dry years than in a normal hydrological year (ε = 0.365 in 2012 vs. ε = 0.032 in 2016) as the structures with high hydraulic conductivity are empty, which increases the water flow velocity. The values for the i function also display high discrepancies In turn, the results of the φ (t) function and the gentle slope of the final section of the hydrograph demonstrate the higher water volume stored in the saturated zone, which causes the steady decrease of the discharge through time. The slow emptying of the phreatic zone (α = 0.005) is responsible for the predominance of the base flow component in the spring outflow (V dyn /V t ), which represents 67.5% of the total water volume drained during the recession period. In dry years, this percentage of base flow increases to 93% (2012) or 87% (2015) of the total water discharged by the spring. The prevalence of base flow in the spring is very high, given that 78% of the days in the recession period have registered a baseflow index above 0.9. The results demonstrate the considerable storage capacity of the saturated zone of the Degracias-Sicó karst aquifer, which can be explained by the large dimensions and Vauclusian characteristics of the saturated zone.
The results of k index range from 0.10 (1.2 months) to 0.20 (2.4 months) with an outstanding peak value in 2012, a remarkably dry year (0.54 = 6.5 months). If one excludes the result of 2012, the residence time of the water drained by OAA spring is relatively short, which could indicate a high karstification degree of the epiphreatic and upper phreatic zones, where a connected system of conduits drains the new infiltrated water to the spring. The development of caves, conduits, and voids in this zone of the aquifer, that simultaneously favors water storage and drainage, could be connected to: (1) the significant fluctuation of the water table from wet to dry seasons (very contrasting in the Mediterranean climate), and (2) the enormous thickness of the Sicó karst massif. This hypothesis could explain the effective regulation power of the aquifer and, at the same time, the large volume of mixed flow drained by the spring during the recession period. The k value in 2012 (0.54) reveals the large storage capacity of the saturated zone of the Degracias-Sicó karst aquifer, responsible for the perennial character of the OAA spring even in very dry years.
According to Mangin's classification of karst systems, the results of the index i and k indicate that the study aquifer is in the domain of complex karst systems, largely extended and made up of several subsystems.

Water temperature and electrical conductivity analysis
The combined analysis of the discharge, water temperature (T) and electrical conductivity (EC) at OAA spring during 3 years, under very different hydrometeorological characteristics, has provided deeper insight and advances on our understanding of the hydrodynamics and structural organization of drainage  (Mangin 1975(Mangin , 1984  patterns inside this aquifer, which was quite unknown from the hydrogeological perspective. The variation of these parameters is displayed in Fig. 6. The thermal behavior of OAA spring exhibits two distinct patterns: the long-term variations, following the seasonal air temperature oscillations, and the short-term variations (few days) as the result of transient recharge episodes.
Analysis of daily data sets has revealed the significant uniformity of the spring water temperature, ranging from 15.7 to 16.7°C, during the study period. The coefficient of variation is very low (1.7%; Table 2), which reflects the long travel time through the aquifer (where an efficient heat exchange occurs between the rock and water) and deep circulation of the water (not influenced by ambient air temperature). Such low variations in the spring water are consistent with the dominance of diffuse recharge and relatively slow groundwater circulation. The recharge water has time to equilibrate to the aquifer temperature, especially in deeper areas of the unsaturated zone.
The thermal regularity of spring water is only interrupted by intense and abundant rainfall events, which means that water traverses the aquifer before reaching thermal equilibrium. The water infiltrates at a different temperature than the aquifer and quickly moves through the system, conserving some characteristics of the input temperature signal. After a heavy rainfall event at during autumn-winter season, the temperature of water drained by OAA spring decreased up to 0.8°C (Fig. 6). This thermal variation is particularly pronounced and rapid after a dry period, since the conduit system was empty. Thus, the magnitude of the deviation of water spring temperature depends on the importance of the recharge event and the previous hydraulic conditions of the system. These short variations indicate the existence of quick circulation through a functional network of conduits and concentrated recharge. The strong negative correlation between the spring discharge and water temperature (r 2 = -0.937, at significance level to 0.01; Table 3) prove that this new water reaches the spring still in the rising phase of the hydrograph.
The electrical conductivity of the spring flow shows similar bimodal behavior to the temperature at seasonal and recharge event scale. The water drained by the spring is significantly mineralized (mean value is 533 μS/cm and maximum values in every year range from 630 to 652 μS/cm). These values of electrical conductivity are very similar to those occasionally registered during the study period at all temporary and perennial springs at the west limit of the aquifer. More than 72% of the daily data present mean values above 500 μS/cm, suggesting a substantial time of residence of the water within the aquifer and, hence, confirms the large storage capacity of the saturated zone and the dominance of base flow in the spring discharge as previously pointed out. The oscillations of electrical conductivity values are not high for a karst system, as the mean coefficient of variation (CV) is 9.8% (Table 2). However, the CV in dry years is much higher (16.2% in 2011/2012), which indicates the existence of karstification and a functional drainage network that quickly drains the infiltrated water to the outlet. López-Chicano et al. (2001) and Liñan et al. (2009) presented CV values of 14 and 19%, respectively, for electrical conductivity data series at springs that drain very karstified systems in the south of Spain.
The higher values of electrical conductivity were registered subsequent to abundant recharge periods and have occurred from autumn to spring season, regardless of the hydrometeorological characteristics of the year. The highest peaks were recorded during the rainiest temporal sequence of the study period (March and April 2013) as the result of significant piston flow mechanisms. The pressure pulse quickly moves through the conduit system, and displaces the water stored in the saturated zone, probably at greater depths with longer residence time than those which usually flow through the saturated zone to the spring. After this initial peak of highly mineralized water, the intense rainfall events and, consequently, the heavy hydraulic charge (due to the rising hydraulic head) continue to activate the flow of more mineralized waters stored in the saturated zone.
The high amounts of rainfall and infiltration at the recharge area, the high thickness of unsaturated zone (around 300 m) and the Vauclusian structure of the aquifer sustain the understanding that the Degracias-Sicó karst aquifer has a very important component of piston flow in the flow regime. The positive correlation between discharge and electrical conductivity (r 2 = 0.227, at significance level to 0.01) shows that part of the rise of spring flow occurs due to the arrival of these highly mineralized waters.
Afterward, the piston flow occurs as a sharp and quick fall of electrical conductivity, which represents the arrival of recently infiltrated water toward the spring. The infiltration water flows vertically through the thick unsaturated zone along the preferential pathways, which increase the flow velocities and inhibit a good mixing process (quick flow regime). Therefore, the spring water is easily diluted during the high flow period and the new water drained by OAA spring (dilution waters) is a mixture of mineralized groundwater and freshwater mainly infiltrated through the concentrated recharge points at the western and south sectors of the aquifer's recharge area.
The progressive rising of electrical conductivity at the hydrograph falling limb phase and the maintenance of those relatively high values (500-600 μS/cm) after an abundant recharge period is linked to a component of slow and diffuse flow from the infiltration surface with longer transit times.  This water, more mineralized than the water from quick flow, corresponds to delayed flow (identified on time series analysis and recession curve analysis), which can feed the spring during several weeks or even a few months. The slow recovery of electrical conductivity and water temperature to the values observed before the recharge event confirms the presence of this type of flow in the spring during the recession period. The combined analysis of the temporal evolution of OAA spring discharge, electrical conductivity and the water temperature has demonstrated the complexity and namely the duality of the hydrogeological functioning of the Degracias-Sicó karst system, including both conduit and diffuse flow behaviors. On the one hand, (1) the coherent but delayed seasonal periodicity of spring groundwater temperature with respect to environmental temperature, (2) the general low variations of thermal and mineral parameters of the spring discharge, and (3) the slow recovery of spring water temperature and electrical conductivity after a peak discharge event, reveal that recharge water moves relatively slowly through the aquifer. This permits an effective heat exchange along the flow paths and, thus, demonstrates the balance of the aquifer characteristics, and this can be interpreted as the result of the lower development of functional karstification in the catchment area drained by this spring and its diffuse flow regime. The relatively stable and high values of both parameters during several months (particularly at low-water conditions) also indicate the considerable storage capacity of the aquifer and the dominance of base flow in the spring.
In contrast, after a rainfall event, the quick and sharp fall of the temperature and electrical conductivity at the OAA spring outflow that occurs immediately after the rising of spring discharge is the result of the arrival of newly infiltrated water from surface recharge points where concentrated infiltration occurs. The short transit time of this new water, especially in dry years, reveals fast circulation throughout the karst system, which proves the existence of a functional karst drainage network. Then, it can infer the effective karstification of the system, which creates flow heterogeneity, increases hydraulic conductivity contrasts, originates different residence times and, consequently, justifies the diverse water-types drained by the spring.

Conclusions
The methodological perspective combining the analysis of the hydrograph and the physical characteristics of OAA spring water has proven to be very useful, reliable and practical to investigate the hydrogeological functioning of a large aquifer and to obtain valuable insights about the internal structure of the karst system.
The results indicate that Degracias-Sicó karst aquifer exhibits a bimodal hydrodynamic functioning. The dominance of a slow flow component which feeds the spring all year, mostly during the depletion period, is associated with the large storage capacity of the aquifer in the phreatic zone and, probably, to the slow and diffuse groundwater circulation through the unsaturated to saturated zone (delayed flow). In turn, the presence of a considerable quick flow component in the spring as the response of intense and abundant rainfall events is associated with the concentrated recharge water circulating through a connected conduit network, which proves the developed karstification of some sectors of this large and complex karst system.
The application of the aforementioned indirect methods has unequivocally displayed the duality of the karst aquifer in terms of: (1) recharge (mainly diffuse but also concentrated into the boreholes); (2) underground circulation (slow and diffuse along the fractures network vs. rapid and concentrated through a structured conduit system); storage (voids, large conduits in the saturated zone vs. rock matrix due to hydraulic gradient inversion during high flow periods). Furthermore, this joint analysis of carefully selected methods has confirmed that statistical components extracted from a discharge signal of a karst spring can provide meaningful information about hydrogeological functioning and physical characteristics of karst systems. The advances presented here are particularly important given the unknown character of hydrodynamic functioning and not-well-known physical features of the study aquifer. Finally, this analysis also highlights the importance of applying several methods in karst studies, as together they permit focus on different areas and aspects of the karst system, particularly the lesser known aspects, where physical data are very difficult or impossible to acquire.
The present study constitutes the first quantitative analysis of the hydrogeomorphological characteristics of the Sicó karst massif through the application of well-defined, rigorous and tested scientific methodologies. Nonetheless, additional effort, particularly hydrochemical analysis and tracer tests, is desirable to get in-depth scientific knowledge and consequently to define specific policies for an effective management of the water resources in this region. you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.