Defining the hydrogeological behavior of karst springs through an integrated analysis: a case study in the Berici Mountains area (Vicenza, NE Italy)

Knowledge of the hydraulic and geological properties of karst systems is particularly valuable to hydrogeologists because these systems represent an important source of potable water in many countries. However, the high heterogeneity that characterizes karst systems complicates the definition of karst hydrogeological properties, and their estimation involves complex and expensive techniques. In this study, a workflow for karst spring characterization was used to analyze two springs, Nanto spring and Mossano spring, located in the Berici Mountains (NE Italy). Based on the data derived from 4 years of continuous hourly monitoring of discharge, water temperature and specific electrical conductivity, a hydrogeological conceptual model for the monitored springs was proposed. Flow rate measurements, which combined recession curve, flow duration curve and autocorrelation function techniques, were used to evaluate the spring discharge variability. Changes in spring discharge can be related both to the degree of karstification/permeability and to the size of the karst aquifer. Moreover, combining monitored parameters and rainfall—analyzed by the cross-correlation function and VESPA (Vulnerability Estimator for Spring Protection Areas) index approach—permitted assessment of the spring response to recharge and the behavior of the drainage system. Although the responses to the recharge events were quite similar, the two springs showed some differences in terms of the degree of karstification. In fact, Mossano spring showed a more developed karst system than Nanto spring. Three systems (two karsts and one matrix/fractured) are outlined for Mossano spring, while two systems (one karst and one matrix/fractured) are outlined for Nanto spring.


Introduction
The analysis of karst springs behavior is a widespread activity in hydrogeology because karst aquifers are a primary source of human water supply in many countries (Ford and Williams 2007). The estimation of available water and the determination of water quality and vulnerability will play a growing strategic role in the future because of the increasing demand for water for drinking, irrigation, and industrial purposes (Assaad and Jordan 1994;Bakalowicz 2005;Crawford and Ulmer 1994;Field 1992;Goldscheider 2012;Kresic and Stevanović 2010). The quantitative characterization of karst systems requires the estimation of realistic hydraulic and geometric parameter values for both the rock matrix and the conduit network (Kovács and Perrochet 2008). Classical geological and hydrogeological surveys of karst (i.e., borehole tests, tracer experiments, speleological and geophysical observations) are not only expensive and sometimes inapplicable but also provide limited information due to the great heterogeneity of karst systems (Kovács et al. 2005).
The monitoring of spring discharge and water physicalchemical parameters (e.g., temperature and electrical conductivity), accompanied by analysis of the resulting hydrographs and chemographs, represents a classical hydrogeological prospection technique that is very helpful when detailed Electronic supplementary material The online version of this article (https://doi.org/10.1007/s10040-020-02122-0) contains supplementary material, which is available to authorized users. analyses such as tracer tests and speleological investigations, are not applicable. This approach allows for assessment of the basic hydrogeological characteristics of a karst aquifer that can be used for planning its utilization and protection. Necessary considerations regarding this method are as follows: (1) the spring hydrograph provides an integrated representation of the network of fractures and conduits transferring groundwater from the recharge to the outflow area and (2) the quality and quantity of groundwater are of diagnostic importance for understanding the functioning of the system (Dreiss 1982(Dreiss , 1989Ford and Williams 2007;Knisel 1972).
Several authors have focused their attention on the analysis of recession curves. Considering a recharge event producing an increase in spring discharge, the recession curve describes the part of the hydrograph showing the progressive decrease in discharge from the peak of the discharge to the pre-storm outflow (Ford and Williams 2007). Starting with the works of Boussinesq (1877) and Maillet (1905), recession curve analysis (RCA) has been used to infer the geometrical or hydraulic properties of aquifers supplying springs (Bonacci 1993;Brutsaert and Nieber 1977;Covington et al. 2009;Dewandel et al. 2003;Fiorillo 2014;Kresic and Stevanović 2010;Padilla et al. 1994;Szilagyi et al. 1998;Troch et al. 1993) or to identify the flow-type induced by different degrees of aquifer karstification (Baedke and Krothe 2001;Birk and Hergarten 2010;Kovács and Perrochet 2008;Kovács et al. 2005;Malík and Vojtková 2012).
Another approach for hydrograph analysis, derived from hydrological practices, is the use of flow duration curves (FDCs). The FDC is a cumulative frequency curve showing the percent of time during which the discharge is equal to or exceeds a specified threshold in a given period. This approach has been used to assess the range and variability of spring discharge (Hartmann et al. 2013;Kovačič 2010;Malík 2015).
Drainage through a karst system can be represented as an impulse function that transforms the recharge input pulse into a flow rate variation at a spring recorded by a hydrograph. The analysis of this function permits the identification of the hydrogeological features of an aquifer and its degree of karstification (Ford and Williams 2007). Time series analysis is a signal processing method employed to explore impulse functions and to improve the understanding of the hydrogeological behavior of a karst system. Mangin (1984) was the first to apply this approach in karst hydrogeology, using autocorrelation and spectral density functions to analyze the spring discharge. Subsequently, different authors have extended the application of these methodologies to other parameters (e.g., precipitation, water temperature, electrical conductivity) by the introduction of other functions such as cross-correlation (Adji and Bahtiar 2016;Bailly-Comte et al. 2011;Fiorillo and Doglioni 2010;Kovačič 2010;Ladouche et al. 2014;Larocque et al. 1998;Mayaud et al. 2014;Zhang et al. 2013). The most widely used methodology consists of the concomitant utilization of autocorrelation functions (ACFs) and crosscorrelation functions (CCFs); in particular, an ACF examines how a value depends on preceding values over a period of time (Larocque et al. 1998;Mayaud et al. 2014). The use of ACFs allows one to quantify the memory effect of a system (Mangin 1984;Panagopoulos and Lambrakis 2006), while in contrast, CCFs provide an overview of the interrelationship between the input and output series (Larocque et al. 1998;Panagopoulos and Lambrakis 2006). In karst hydrogeology, a CCF is usually employed to evaluate the correlation between rainfall and spring discharge, highlighting the response of a system to the recharge impulse (Chung et al. 2018;Fabbri et al. 2013;Fiorillo and Doglioni 2010;Liu et al. 2011;Tagne and Dowling 2018;Zhang et al. 2013).
The Vulnerability Estimator for Spring Protection Areas or VESPA index (Galleani et al. 2011) is a promising method used to assess the degree of vulnerability of a karst spring, but it can also be used to classify general hydrogeological behavior through the analysis of spring responses to recharge impulses. The prevailing behavior is identified considering two parameters: spring discharge and electrical conductivity (Banzato et al. 2017). Moreover, an advantage of this method is the ability to combine these two parameters with water temperature to assess spring vulnerability.
In this study, the RCA, FDC, ACF and CCF techniques were used to define the conceptual hydrogeological model of two springs located in the Berici Mountains (northeastern Italy). In addition, the parameters defined by the VESPA index were used to assess the type of flow and vulnerability of springs. Moreover, particular attention was paid to the meaning of the CCF analysis between spring discharge and electrical conductivity and between spring discharge and water temperature as an estimate of the minimum transit time inside the conduit network feeding the springs. These techniques enable use of only the monitoring parameters together with rainfall events to obtain qualitative information about the hydrogeological features characterizing the analyzed springs.
The objectives of this study are twofold: (1) outline a workflow for karst spring characterization and (2) contribute to highlighting the main features of groundwater circulation in a peculiar karst system. In fact, the Berici Mountains represent a peculiar karst system because they appear as a calcareous plateau, surrounded by the Venetian Plain and isolated from the nearby Prealps.

Geological and hydrogeological settings
The Berici Mountains extend across the central part of the Veneto Region to the southwest of the city of Vicenza (Veneto,NE Italy;Fig. 1a). Together with the Lessini Mountains and the Euganei Hills, the Berici Mountains constitute a relatively undeformed foreland structural high known as the Lessini-Berici-Euganei block (LBE; Pola et al. 2014). Toward the east, the triangular-shaped LBE block is separated from the Eastern Southern Alps and their deformed foreland by the high-angle NE-dipping Schio-Vicenza fault (NW-SE trend; Fig. 1a). Within the LBE block, the Berici Mountains are separated from the Euganei Hills to the south by the high-angle Riviera Berica fault (NNE-SSW trend). These faults are buried beneath the Quaternary alluvial deposits (Benvenuti and Norinelli 1967;Mietto 1988a;Pola et al. 2014).
The Berici Mountains are an isolated plateau covering almost 200 km 2 (Fig. 1b) with a maximum elevation of approximately 450 m above sea level (asl) (Monte Alto, elevation 444 m asl). Two long and deep valleys, Val Liona and Valle di Fimon (Fig. 1b), extend to the core of the complex, which can then be divided into two distinct sectors, the eastern and western sectors (Mietto 1988a;Mietto and Sauro 2000). The stratigraphic sequence is characterized by Upper Cretaceous to Miocene sedimentary rocks with subhorizontal or gently dipping beds and by basaltic rocks related to the Paleogene Lessini-Berici-Euganean volcanic cycle (Mietto 1988a).
The study area is located in the eastern sector of the Berici Mountains (Fig. 1b). The oldest outcropping rocks are represented by upper Cretaceous marly limestone of the Scaglia Rossa Formation (SR), with a thickness of some tens of meters, followed by middle Eocene marly arenaceous limestones (ML, maximum thickness of 35 m) and Nummulitic limestones (NL,   In legend, information about the permeability degree of each unit are presented (Fabiani 1911;ARPAV 2007). The coordinates are in UTM zone 32 N system using the WGS84 datum (expressed in km) that can reach a maximum thickness of approximately 200 m. The sedimentary sequence is closed by the Oligocene limestones of the Castelgomberto Formation (CFm), with a thickness ranging from 80 to 200 m. The whole succession was affected by a basaltic explosive volcanism (Oligocene), with the formation of several necks, volcanoclastic rocks and tuffites of variable thickness (BV; Fig. 2 and Table 1; Bassi et al. 2000;Mietto 1988a). Karst phenomena played a significant role in the geomorphological evolution of the area, leading to the formation of both surficial and underground karst features. This phenomenon is particularly true for the high-karstified limestones of CFm, where the dolines are the most representative large-scale karst features (Castiglioni 1996;Mietto 1988b;Mietto and Sauro 2000;Sauro 2002Sauro , 2005. Considering the hypogean features, approximately 630 caves are present in the Berici Mountains, many of which are located in the CFm (Fig.  1). The spatial extension of these caves ranges between a few meters and approximately 1 km, while the maximum difference in elevation is approximately 100 m ( Fig. 2 Geologic map of the study area. Lithologies were obtained from the lithological map of the Veneto Region (1:250,000; Antonelli et al. 1990), while the stratigraphic contacts between the different units were detailed during field work. From the bottom: Scaglia Rossa Formation (SR), Pelagic limestone containing planktonic foraminifera; Marlyarenaceous Limestone (ML), alternating marly limestone, calcareous marls, and calcarenites interbedded with tuffs; Nummulitic Limestone (NL) massive limestone rich in Nummulites; Priabona Formation (PFm), marly limestones with rich fossil content; Castelgomberto Formation (CFm), massive or irregularly bedded limestones and irregularly bedded calcarenites with miliolids and coralline red algae; basaltic volcanits (BV), basaltic explosive volcanism represented by necks, volcanoclastic rocks and tuffites (Bassi et al. 2000). Regarding settlements: Nt Nanto, Ms Mossano. In the top-left box is a rose diagram of the orientations of the maximum diameter directions of dolines (DMAX; Torresan 2016). The coordinates are in UTM-zone 32-N system using the WGS84 datum (expressed in m) because of their degree of karstification, while SR and ML can be considered aquitards due to their greater terrigenous composition. In addition, due to the low permeability of the volcanic products, BV can be considered a permeability barrier for the groundwater flow affecting the karst systems (ARPAV 2007;Fabiani 1911).
As typically occurs in karst environments, dolines represent a preferential path for the infiltration of meteoric water, even though their bottom can be infilled with fine sediments. As a consequence, the Berici Mountains landscape is characterized by the absence of any hydrographic network, and its recharge is exclusively of autogenic origin (Ford and Williams 2007). Diffuse autogenic recharge is channelized in a welldeveloped underground karst system that allows the circulation of infiltration waters through different-sized conduits, fractures, etc. Subsequently, water discharges from several springs located at the base of the Berici plateau and preferentially along the contact between permeable rock units (NL, PFm or CFm) and low-permeability layers (ML or SR;Fabiani 1911). This system generates a large number of springs, whereby approximately 166 springs are present in the Berici Mountains (Marcolongo 2005). Among these springs, 43 springs show a permanent regime, while the remaining 123 springs are periodic (periods without flow during the long dry season) or episodic (periods of activity in response to only main rainfall events). Most of the permanent springs have an average discharge that is generally lower than 2 L/s, while the maximum average discharge is approximately 20 L/s (Marcolongo 2005).
The monitored springs, called Nanto and Mossano, are two perennial springs located in the eastern part of the Berici Mountains in municipalities of the same names (Fig. 1b). The decision to monitor these two springs, among all the springs located in the Berici Mountains, is related to some pre-established requirements: (1) perennial regime; (2) accessibility; (3) presence of an intake capture structure; and (4) lack of water exploitation.
Nanto spring (Fig. 3a) discharges at the base of the slope of the Berici Mountains at an elevation of approximately 25 m asl The spring emerges at the contact between NL and the underlying ML (Figs. 2 and 3b).
Mossano spring (Fig. 3c), located at an elevation of approximately 90 m asl, also discharges at the contact between NL and ML (Figs. 2 and 3d). Information about the monitored parameters for both springs are shown in the 'Results and discussion' section.
Dolines in the karst area surrounding the two springs were recognized by the analysis of satellite images (Landsat) and a digital elevation model (DEM) in a geographic information system (GIS) environment (Fig. 2). The alignment of several dolines leads to the assumption that some tectonic structures, i.e., fractures and faults, exist. The interaction between tectonics and karstification was confirmed upon considering the orientations of the maximum diameter directions of dolines (DMAX; Bondesan et al. 1992). As shown in Fig. 2, the DMAXs are preferentially WNW-oriented according to the direction of the identified structures (Torresan 2016). In addition, a structural field analysis performed near Mossano spring allowed these elements to be characterized as extensional faults related to Paleogene extensional tectonics (Zampieri 1995).

Monitoring data
Nanto and Mossano springs were continuously monitored for almost 4 years, from January 2015 to June 2018, by two highaccuracy multiprobe sensors (CTD-Diver). The recorded parameters were (1) water-level height above the diver (H in cm; range 0-1000 cm, accuracy ±0.5 cm and resolution 0.2 cm), (2) temperature (T in°C; range from −20 to +80°C, accuracy ±0.1°C and resolution 0.01°C), and (3) specific electrical conductivity at 25°C (EC in mS/cm; range 0-30 mS/cm, accuracy ±1% of reading and resolution ±0.1% of reading). An acquisition time interval of 1 h was chosen to register the effect of flood events during the hydrogeological year (Galleani et al. 2011). The water level measurements were converted to discharge (Q) values by two experimental level-discharge relations obtained through 57 flow-rate and level measurements collected from January 2015 to June 2018 (see the electronic supplementary material, ESM).
The time series of both springs were subject to several interruptions due to data logger malfunctions. The main monitoring interruptions for Nanto spring occurred from September 2015 to January 2016, from October 2016 to January 2017 (only for EC), and from February 2017 to April 2017. Mossano spring experienced only one long interruption, from September 2015 to January 2016. Moreover, both springs exhibited short interruptions during the monitoring period. It is worth mentioning that the CTD-Diver used in Mossano spring was located in an external storage tank (Fig.  3c), in which the water temperature was affected by atmospheric agents and by daily and seasonal air temperature variations. Conversely, the measurement of the water level was not directly influenced by rainfall events in the spring area. For this reason, a new data logger for Mossano spring was added on January 2018 with the aim of measuring only water temperature. This data logger was positioned in the collection chamber connected to the external storage tank (Fig. 3c), where it was impossible to conduct water level measurements. Data collected by the Barbarano Vicentino meteorological station were used to evaluate the local precipitation and temperature and to perform the time series analysis. This station is part of a public monitoring network covering the Veneto Region, with 167 meteorological stations managed by the Regional Agency for Environmental Prevention and Protection (ARPAV). The Barbarano Vicentino meteorological station is located at an elevation of 16 m asl, 2.9 km from Mossano spring and 4.2 km from Nanto spring (Fig. 1b), representing the nearest station to the monitored springs. Although the infiltration area of the two karst systems is located at an elevation of approximately 300-400 m asl, there are no meteorological stations representative of these elevations close to the Berici Mountains.

Time series analysis and vulnerability
The workflow used to characterize the two monitored springs, starting from the monitored parameters, is schematized in Fig. 4. The analysis of the variations in spring discharge and the definition of the behavior of the springs were used to delineate the conceptual hydrogeological model.
The discharge variations recorded in karstic springs are mainly related to the degree of recharge and karstification of the system. For this reason, a detailed analysis of spring discharge must be performed to assess some of the main features characterizing a given type of aquifer. In the proposed workflow, three methods were chosen to analyze the flow rate variability: (1) recession curves analysis (RCA), (2) flow duration curve (FDC) analysis, and (3) autocorrelation function (ACF) analysis.
The spring response in relation to rainfall events indicates the principal characteristics of the monitored karst system. To evaluate the spring behavior, the variations in the other monitored parameters (Q, T and EC) and their relationships to recharge were considered. Two types of analyses were performed: cross-correlation functions (CCFs) and VESPA index analysis.
Recession curves analysis allows for the study of variations in spring discharge as a consequence of recharge events. The analysis of the recession curves focused on the recession coefficient (α), which is qualitatively related to the degree of aquifer karstification (Baedke and Krothe 2001;Birk and Hergarten 2010;Dewandel et al. 2003;Fiorillo 2014;Kovács and Perrochet 2008;Kovács et al. 2005;Malík and Vojtková 2012;Tallaksen 1995) and, in part, to the size of the karst system (Bonacci 1993). The recession coefficient was used by Maillet (1905) to describe flow recession using a simple exponential equation (Bonacci 1993;Ford and Williams 2007): where Q t is the discharge at time t, Q 0 is the initial discharge at time t 0 and α is the recession coefficient. Boussinesq (1903Boussinesq ( , 1904 express the recession by a quadratic form: Dewandel et al. (2003) pointed out that Eq.
(1) provides an approximate analytical solution for the diffusion equation that describes flow in a porous medium, whereas Boussinesq's quadratic equation provides an exact analytical solution (Ford and Williams 2007).
However, the mathematical properties of Eq. (1), in which the derivative is exponential and α is independent of Q 0 , are more convenient than the quadratic solution in Eq. (2), explaining the widest use of the Maillet formula (Dewandel et al. 2003).
Although Eq. (1) is useful for describing the baseflow recession of a karst system (Kovács and Perrochet 2008;Kovács et al. 2005; Malík and Vojtková 2012), many authors have proposed describing the cumulative effects of several reservoirs in order to interpret the entire recession curve (Amit et al. 2002;Dewandel et al. 2003;Ford and Williams 2007;Kovács and Perrochet 2008). As a matter of fact, the recession curve can be subdivided into multiple segments, describing different reservoirs with different hydrogeological properties (i.e., permeability, porosity, storativity), all contributing to spring discharge. For each segment, there is a corresponding α value obtained by the exponential Eq. (1). According to Taylor and Greene (2008), a difference in α exceeding one order of magnitude can be associated with different types of flow regimes and is used to describe spring behavior during the recession period. High values of α indicate fast drainage (conduit dominated flow) related to well-developed karst conduits and, thus, to a high degree of karstification. Conversely, low values of α are ascribable to diffuse flow into smaller fractures and pores (diffuse dominated flow), indicating a lower degree of  Fig. 4 The workflow adopted in this work to achieve the final definition of a hydrogeological conceptual model. Flow rate data are used to evaluate the spring discharge variations through recession curve analysis (RCA), flow duration curve (FDC) and auto-correlation function (ACF). The monitoring parameters (Q, T and EC) and the rainfall events (R) are used to evaluate the spring behavior by the cross-correlation function (CCF) and VESPA index approach. For each approach is defined the main parameter evaluated and the types of information that can be achieved karstification (Bagheri et al. 2016;Bonacci 1993;Ford and Williams 2007;Kovács and Perrochet 2008). Flow duration curves analysis is a measure of the range and variability of spring discharge. FDC analysis represents the percentage of time during which the flow rates of a given spring exceed a specified threshold (Malík 2015). The discharges were ordered and plotted from the highest to the lowest values, without considering their chronological sequence. The obtained graph shows the discharge on the ordinate axis and the percentage of the monitoring time that exceeded the specified flow-rate threshold on the abscissa axis. FDC can be used to estimate principal karst system features by analyzing the slope of the curve. In fact, the presence of high slopes suggests well-developed karst conduits (conduit dominant flow), whereas low slopes indicate systems with a high storage capability and a lower degree of karstification (diffuse dominant flow; Malík 2015).

Methods Parameters Information
Moreover, the use of the autocorrelation function permits the definition of the memory of a karst system in terms of discharge variations. The ACF can be defined as a normalized measure r(k) of the linear dependence among successive values of the data series (discharge values in this case), and it is as described (Zhang et al. 2013): in which: where k is the time lag, n is the length of the time series, t is the time, x is a single event measurement, and x is the mean of the event measurements. The related autocorrelogram is useful for evaluating the memory of a system, which is estimated through the decorrelation lag time. This lag time is defined as the time needed by the autocorrelation function to reach a determined value, usually between 0.1 and 0.2 (Mangin 1984;Liu et al. 2011;Panagopoulos and Lambrakis 2006;Tagne and Dowling 2018;Zhang et al. 2013). The autocorrelation values of an uncorrelated time series rapidly decrease (Chung et al. 2018;Liu et al. 2011), which can be interpreted as a system with short memory. Conversely, a high system memory results in a gentle decrease in the slope of the autocorrelation function; therefore, through autocorrelogram analysis, some consideration related to the degree of karstification and the size of the karst aquifer can be made (Liu et al. 2011;Mangin 1984;Padilla and Pulido-Bosch 1995).
To evaluate the springs' responses over time with respect to recharge, CCFs were carried out with the monitored spring parameters (Q, T, EC) and rainfall (R). Cross-correlation is widely used in hydrogeology to analyze linear relationships between input and output signals (Chung et al. 2018;Fabbri et al. 2013;Fiorillo and Doglioni 2010;Liu et al. 2011;Tagne and Dowling 2018;Zhang et al. 2013).
This statistical analysis was performed using R code (R Core Team 2018). Before the application of the CCFs, the data series were filtered because the parameters Q, T and EC are usually autocorrelated. The filtering process was performed with the use of a simple filter and defined as the differences of one lag for every parameter. This approach corresponds to determining the difference between each value and its previous value, removing a trend and making stationary a nonstationary time series (R Core Time 2018). CCF represents the correlations between two different data series with respect to lag times. The cross-covariance function, where u is the input and y is the output, is given by Eqs. (6) and (7): where u and y are the averages of the input and output parameters, respectively. CCF is defined as: Equation (7) was used primarily to analyze the crosscorrelation between R (input) and Q (output), to define the system response and to qualitatively evaluate the travel time of the meteoric water to cross the vadose zone. When there was a correlation to a certain lag, this was interpreted as the delay (expressed in hours) of the spring response to a recharge event. Furthermore, cross-correlations between spring parameters, i.e., Q vs. T and Q vs. EC, were performed to assess the spring behavior. Q vs. EC allows the evaluation of the minimum travel time of the neo-infiltrated water to transit through the conduits system and reach the point of emergence at the springs. Consequently, information about the size of the karst system can be obtained. Additionally, T vs. EC was calculated to evaluate whether these two parameters were synchronous. The CCF results were plotted in two-dimensional (2D) graphs, known as cross-correlograms, in which the correlation value (between −1 and +1) is reported on the y-axis, and the lag is shown on the x-axis. In addition, two confidence bands were displayed to graphically define whether or not the two analyzed series are correlated. If two series are correlated, the correlation value must be higher than the confidence bands. These bands, one negative and one positive, were defined as: − 1:96 ffiffi n p and 1:96 ffiffi n p , respectively (R Core Team 2018).
Through the use of more than 1 year of continuous hourly monitoring data, the type of spring behavior was classified using parameters of the VESPA index (Galleani et al. 2011) based on the Q, EC and T values of the spring. The VESPA index V was defined by the following relationship (Banzato et al. 2017): In Eq. (9), c(ρ) is the correlation factor, ranging from 1 to 0, defined as: where ρ is the correlation coefficient between Q and EC, ranging between −1 and 1, calculated over more than 1 year of continuous hourly data. u(ρ) is the Heavside step function: where α is a scale coefficient ranging from 0 and 1, which is normally assumed to be 0.5, and β and γ are the temperature and discharge factors, respectively, defined as: where T max and T min are the maximum and minimum temperatures, respectively, and Q max , Q min and Q med are the maximum, minimum and average discharge values recorded during the monitoring period, respectively. The use of Eq. (9) allows one to define the spring vulnerability degree, as shown in Table 1. Because Eq. (9) is related to the joint use of Q, T and EC, each factor can increase or decrease the spring vulnerability. No recharge measurements are required for the definition of the VESPA index, and there are no relationships between recharge and discharge. The four degrees of vulnerability shown in Table 1 are the result of experimental applications of the VESPA index and local regulations (Piedmont Region, NW Italy), defining the vulnerability of a karst spring area (Galleani et al. 2011).
By means of ρ and c(ρ), Galleani et al. (2011) proposed three different types of prevailing spring phenomena (Table 2): 1. The type A "replacement" phenomenon represents a welldeveloped discharge system. The replacement effect produces an increase in Q, a variation in T according to the seasonal rainfall temperatures and a decrease in EC. This replacement effect is typical for highly karstified limestone aquifers. 2. The type B "piston" phenomenon represents a moderate drainage system in which the hydrodynamic response can show impulse behavior. This phenomenon is typical for slightly karstified limestone aquifers with significant water storage. The piston effect produces increases in Q, T and EC. 3. The type C "homogenization" phenomenon is characterized by a poorly developed drainage system, where the water response to rainfall is almost absent. The Q trend presents slow and modest fluctuations that are uncorrelated with T or EC.
Due to the differences in the development of the fracture network, the type A drainage systems are generally the most vulnerable, while type C drainage systems are the least vulnerable (Galleani et al. 2011).
The type of drainage system is inferable on the basis of c(ρ) and ρ, as in Table 2. In addition, because the coefficient ρ is defined as the correlation between Q and EC, a graphical comparison between these two parameters was performed to validate/deny the behavior classification and to explain the EC variations that are ascribable to the groundwater residence time in the karst aquifer. For this reason, the graphical comparison between Q and EC was not included in the workflow in Fig. 4 because it represents a supplementary analysis.

Results and discussion
The analyzed monitored period utilized to perform the statistical analyses was from April 2017 to June 2018. This period represents the longer-interval time without significant  (Table 3) and histograms (Fig. 5) of the monitoring data, it is possible to highlight a difference in the distribution of data between Nanto and Mossano springs. In terms of discharge, the main differences are as follows: (1) the Mossano Q distribution is bimodal, while the Nanto Q distribution is unimodal; (2) the maximum value of Nanto is higher than the maximum value of Mossano; and (3) the median and mean are higher in Mossano than in Nanto. Likewise, the EC values of Nanto are bimodal, and the median and mean values in Nanto are higher than those in Mossano. Furthermore, the Mossano T values are lower than the Nanto T values, and both springs show a unimodal histogram distribution of T. Regarding the rainfall data, Table 4 shows a summary of the parameters recorded over the 4 years of monitoring by the Barbarano Vicentino meteorological station.

Spring discharge variations
Maillet's equation was used to calculate the recession coefficient α for different segments of the analyzed discharge  . For this type of analysis, compared to the statistical method, the entire monitored period can be used because RCA is less influenced by data continuity. Two representative recession curves among those analyzed (see ESM), with the related α values, are shown in Fig. 6. The α values calculated for the Nanto and Mossano springs highlight the different karstification degrees of the lithologies forming the aquifer reservoir. The first system, described by α coefficients within 10 −1 order of magnitude, is constituted by well-developed conduits, allowing for the rapid discharge of rainfall infiltration (conduit dominant flow; Taylor and Greene 2008). This system is probably more developed in the limestones of the CFm and NL formations than in the marly limestone of PFm. The second system is characterized by recession coefficients within 10 −3 orders of magnitude, suggesting a slower circulation in fractured rocks (diffuse dominant flow; Taylor and Greene 2008) involving both NL. and PFm. The main difference between Nanto and Mossano is the existence of a recession value of 10 −2 for Mossano spring, suggesting an intermediate or mixed flow system (i.e., a system bearing both karstic and fissured characteristics). This hydrogeological behavior represents a transition between the conduit system and the matrix-fracture   Taylor and Green (2008) system and could be related to a karstification degree of PFm that is possibly greater in the Mossano recharge area than in the Nanto recharge area, as testified by field observations (Torresan 2016). The calculated α values (Fig. 6) are greater than the recession coefficients that are commonly available in the literature-for example, the recession coefficients found by Amit et al. (2002) range from 0.1 to 0.0006, while Bagheri et al. (2016) showed α values ranging between 0.02 and 0.0002; furthermore, in Giacopetti et al. (2017), the recession coefficients ranged between 0.01 and 0.001. The differences in these α coefficients are probably due to both the geological and hydrogeological features of the recharge areas and to the small extension of the karst systems. As a matter of fact, as stated by several authors (Baedke and Krothe 2001;Bonacci 1993;Birk and Hergarten 2010;Fiorillo 2014;Kovács and Perrochet 2008;Kovács et al. 2005), the recession coefficient is inversely proportional to the area of the active aquifer and/or to the recharge area feeding the spring. The type of drainage system suggested by the RCA is also confirmed by the FDCs (Fig. 7). The steep slope of the upper end of the FDCs indicates the existence of karst conduits, resulting in a highly variable spring discharge. Conversely, the semiflat slope at the lower end of the curves indicates high water storage due to the low-permeability system (fracture-matrix system). In addition, differences in the behavior between the two springs are highlighted. Nanto spring shows essentially two components, while Mossano spring shows three components. This difference can be related to the existence of an intermediate system (mixed flow, shown in Fig. 6) with a remarkable capability of discharge (9.4-10.4 L/s, exceedance of 10-50%).
ACF analysis further confirms the different behaviors of the two monitored springs. The autocorrelogram analyses (Fig. 8) highlight an initial component characterized by a steeper slope. The ACFs quickly decrease in the first 50-70 h, indicating the existence of well-developed karst conduits that favor a rapid discharge flow. A second component is present only in the Mossano ACF (Fig. 8b) and shows an intermediate slope with a moderately slow reduction of the ACF ranging from approximately 50-900 h. This component can be ascribed to the existence of an intermediate flow system that is also distinguished by RCA (Fig. 6) and FDCs (Fig. 7). Finally, the third component, which reaches the threshold value of 0.2, developed from approximately 50-320 h (almost 13 days from lag 0) in Nanto spring and from approximately 900-1,000 h (approximately 40 days from lag 0) in Mossano spring. This third component confirms the presence of a low-permeability fracture network, which is more extensive in Nanto than in Mossano, that influences the storage capability, delaying the groundwater flow toward the springs.
The lag times required to reach the threshold value allow for some consideration of the size of the karst systems. When the flow rate in Nanto spring reaches the threshold value, Mossano spring has an autocorrelation of approximately 0.5. This difference can be related to the size of the hydrogeological system, which is probably greater in Mossano than in Nanto.

Spring behavior
To evaluate the spring responses, cross-correlation (CCF) analyses were performed considering the period from January 2018 to June 2018 (Figs. 9 and 10). This interval time represents the monitoring period in which the parameters Q, T and EC are representative for both springs and without missing data. The cross-correlation between R (input) and Q (output) allows one to identify the spring response in terms of discharge variations as a consequence of recharge events. Figure 9a clearly shows positive correlations at lags of 4 and 5 h for Nanto spring, meaning a response delay of discharge occurs 4-5 h after a rainfall event. Additionally, for Mossano spring, the significant positive correlations at lags of 5 and 6 h (Fig. 10a) define a response occurring with a delay of 5-6 h after a rainfall event. However, both springs also showed negative correlations (lags of 1 and 2 h for Nanto spring and lags of 3 h for Mossano spring). A possible interpretation of this result is that 1-3 h after rainfall events (increased R), the discharge is decreasing due to its natural trend before the recharge flow perturbation. This behavior is interpreted as the necessary travel time for neo-infiltrated water to cross the vadose percolation zone and reach the saturated zone of the karst aquifer, increasing its hydraulic head and, as consequence, the spring discharge. The CCFs of Q vs. T and Q vs. EC are also useful for analyzing the springs' behavior. Both springs showed only negative significant correlations (Figs. 9b,c and 10b,c), indicating a contrasting response between discharge and both temperature and specific electrical conductivity, traducing in a replacement behavior (type A drainage system). In particular, CCFs between Q and T indicate a maximum negative crosscorrelation at a 1-h lag in Mossano and a 4-h lag in Nanto (Figs. 9b and 10b). Similarly, the CCFs between Q and EC show a maximum negative correlation at a 2-h lag in Mossano and a 4-h lag in Nanto (Figs. 9c and 10c). Of interest are the differences in the response delay between the input (Q) and output (T, EC). The Nanto negative correlations for T and EC have 2-6 h delays in comparison with the discharge (Fig.  9b,c), while they have 1-3 h delays in Mossano (Fig. 10b,c). Furthermore, the CCFs between T and EC showed a very high positive correlation at lags of 0 and 1 h ( Fig. 9d and Fig. 10d), showing that these two parameters behave in the same way. vs T -Nanto Q vs EC -Nanto T vs EC -Nanto Fig. 9 Cross-correlation between the monitored parameters in Nanto spring in the period between January 2018 and June 2018. a Crosscorrelation between R vs Q. b Cross-correlation between Q vs T. c Cross-correlation between Q vs EC. d Cross-correlation between T vs EC. The dashed blue lines represent the confidence bands emergence (e.g., the first arrival time of a tracer test; Ford and Williams 2007). Considering that EC is less subject to external agents compared to T, it is possible to demonstrate the different hydrogeological behaviors of the investigated system. In Nanto, the first negative correlation between Q and EC at lag 0 can be related to the fact that once the freshwater has infiltrated the karst aquifer, it reaches the spring in less than 1 h. All the negative correlations between lags 0 and 7 can be interpreted as the time needed for all of the neo-infiltrated water to reach the spring, with a maximum peak of concentration at 4-h lag. In Mossano spring, the first negative correlation is at a 1-h lag, indicating that after infiltration into the karst reservoir, the freshwater reaches the emergence point after 1 h. The total time required to deplete the neoinfiltrated water is 3 h, with a maximum peak at 2-h lag. Two main differences can be observed between Nanto and Mossano springs: (1) the instantaneous response in Nanto can be related to both the existence of a well-developed karst system and the size of the hydrogeological system (smaller in Nanto than in Mossano), and (2) the longer time needed to deplete the freshwater in Nanto is probably due to a high contribution of a fracture/matrix system rather than a conduit system. Conversely, the shorter time in Mossano may be related to the concomitant influence of three systems (instead of two), namely, the conduit system, the intermediate system and the fracture/matrix system, with the intermediate system developed inside the marly limestones of the PFm. The VESPA index was calculated considering more than 1 year of continuous hourly monitoring data (April 2017-June 2018), without a large number of missing data, as suggested by Banzato et al. (2017). The results are summarized in Table 5. Nanto is classified as a highly vulnerable spring (V = 1.78; Table 1), which exhibits with homogenization behavior (ρ = 0.14; c(ρ) = 0.07; type C in Table 2). Conversely, Mossano is classified as a replacement spring (type A in Table 2) by the coefficient ρ and as a replacement/piston (type A/B) spring by the factor c(ρ), while the VESPA index (V = 13.04) points to a very high degree of vulnerability (   Nanto spring and Mossano spring. Because the spring behavior is defined as the correlation between Q and EC, it is worth observing the related hydro-and chemo-graphs (Fig. 11). Both springs exhibit a comparable discharge trend (Fig. 11a), which is supported by the positive correlation of 0.66; however, Mossano spring shows a greater variability in the flow rate than that of Nanto spring. These discharge variations can explain the higher degree of vulnerability of Mossano spring (Table 5). Regarding the EC, both springs shown decrease in EC, corresponding to the main Q peaks (Fig. 11b), indicating a replacement behavior, as suggested by the CCFs analyses. This observation highlights the existence of well-developed karst conduits that rapidly discharge infiltrated waters and that are poor in ion content, as a consequence of recharge events. In addition to this behavior, Nanto spring show a general increasing trend, which is more prominent after the time 4,000 h (Fig. 11b), that is related to an increase in the ion content of the groundwater due to a lowpermeability system (fractures/matrix system). In particular, the increase in EC is more accentuated after the EC negative peaks for Nanto spring. In terms of spring behavior, the conduit system shows a replacement behavior, while a homogenization behavior can be attributed to the fractures/matrix system. The VESPA index does not allow the separation of these two contributions; however, the results obtained by the analysis of the ρ and c(ρ) factors for Nanto spring emphasized that the homogenization drainage system prevails over the replacement system, confirming that the fractures/matrix system is more extensive than the karst system. In contrast, excluding the EC negative peaks corresponding to the recharge events, Mossano spring shows a generally constant EC trend (Fig. 11b). In addition, after the decrease in EC corresponding to the Q peaks, the EC recovery is gradual. The difference between Nanto and Mossano springs can be explained by the presence of an intermediate system in Mossano, between the conduits system and the fractures/ matrix system, producing a progressive transition from high flow to baseflow. The progressive recovery of EC values in Mossano spring defines three permeability systems that are responsible for the general replacement behavior individuated by the VESPA index analysis.

Hydrogeological conceptual model
According to the proposed workflow (Fig. 4) Fig. 11 Comparison between the two analyzed spring parameters in the period between April 2017 and June 2018: a comparison between Nanto discharge and Mossano discharge. b Comparison between Nanto specific electrical conductivity and Mossano specific electrical conductivity conceptual model for both springs can be delineated. The analysis of the variations in spring discharge, accomplished using the RCA, ACF and FDC techniques, allow one to obtain qualitative information about the permeability of the karst system, related to the degree of karstification, and to the size of the karst system. The study of the spring behavior, accomplished using the CCFs and VESPA index analyses combined with monitored parameters (Q, T, EC) and the recharge (R), allows information to be obtained about: (1) the spring response, which can be related to the neo-infiltrated water travel time across both the vadose zone and the saturated zone, (2) the size of the karst systems and (3) the degree of vulnerability of the karst springs. Conceptually, it is possible to consider the two analyzed karst systems as a multi-reservoir system to explain their hydrogeological structure. In the case of Nanto spring, two reservoirs with different karstification degrees can be individuated. In fact, RCA, FDC and ACF analyses delineated the existence of both well-developed karst conduits and a fractures/matrix system of lower permeability. In addition, the fractures/matrix system is more developed in size than the conduit system, as indicated by the FDC and ACF results. After a rainfall event, the neo-infiltrated water requires 4-5 h to move across the unsaturated zone and reach the karst aquifer, as shown by the CCF between Q vs. R. Subsequently, freshwater reaches the point of emergence of the spring through the well-developed karst conduits in less than an hour (Q vs. EC). However, the contribution of karst conduits to the depletion of neo-infiltrated water is less than that of the fractures/matrix system, shown by the CCF between Q vs. EC. The abrupt increases in the EC after the EC negative peaks (Fig. 11b), in conjunction with the recharge events, are attributable to the net transition between a conduitdominant flow to a diffuse-dominant flow.
The results obtained through the use of RCA, FDC and ACF analyses allowed for the identification of the existence of three reservoirs affecting Mossano spring. In particular, a karst conduit system and a fracture/matrix system, separated by an intermediate system, were individuated. The intermediate system is relevant in the existence of the Mossano spring. In fact, the intermediate system is active for a longer time compared to the other two systems, as indicated by the ACF and FDC analyses. The intermediate system, constituted by both karstic and fissured features, is more developed in size compared to the other two systems. During a recharge event, the infiltrated freshwater reaches the karst aquifer after 5-6 h, as indicated by the CCF between Q vs. R. Due to the high permeability of the karst conduits, the travel time of the neoinfiltrated water across the karst reservoir is shorter (about 1 h) and the water is depleted after approximately 3 h (the CCF between Q vs. EC). The existence of the intermediate system allows for the gradual recovery of the baseflow stage, as suggested by the progressive increase in EC after the negative peaks, related to infiltration processes (Fig. 11b).
Comparing Nanto spring and Mossano spring, some differences appear. In terms of the degree of karstification, the Mossano system is more karstified than the Nanto system. In fact, in the Nanto system, the fracture/matrix system prevails over the conduit system, while, in Mossano, an intermediate system exists. This conclusion is further emphasized by the results of the RCA, FDC and ACF analyses. The intermediate system affecting Mossano spring is probably due to the PFm that may be more permeable in the Mossano karst area.
Regarding the size of the karst domain, ACF analysis shows that the Mossano karst system is greater than the Nanto system. Additionally, the CCF between Q vs. EC confirms the different size by comparing the time of the first arrival of freshwater at the springs. However, the transit time needed for the water to infiltrate the epikarst is comparable (4-6 h). This last aspect is amenable to the same features of the recharge area, constituted by an alteration of the CFm with the formation of soil, mainly used for agricultural purposes.
Finally, other differences are highlighted by the VESPA index that identify Nanto as a highly vulnerable spring with a type C drainage system, while Mossano as a very highly vulnerable spring characterized by a type A drainage system.

Conclusions
In this study, a workflow for the hydrogeological characterization of karst springs was proposed. This workflow was applied to the study of two springs located in the Berici Mountains (Vicenza, NE Italy), with the aim of increasing the knowledge of groundwater circulation in a peculiar karst system.
A multi-approach for the analysis of spring discharge variations, using recession curve analyses, flow duration curves and auto-correlation function techniques, was employed to evaluated the karst degree and the size of the system. Moreover, by combining the monitoring parameters (Q, T, EC) with recharge events (R) through the joint use of cross-correlation functions and the VESPA index approach, it was possible to obtain qualitative information of: (1) the spring response to recharge events, (2) the travel time of meteoric water into both the vadose zone and the conduits network (3) the size/permeability of the karst systems, (4) the spring vulnerability and (5) the drainage system behavior. The information obtained by all the performed analyses allows one to define a hydrogeological conceptual model for each spring that can be used to explain the main features that characterize each karst system. In fact, the previously mentioned quantitative/qualitative information comprises peculiar points in the preliminarily studies of a karst spring, finalized to its potentiality assessment in the water supply.
The two monitored springs showed the same behavior as a consequence of the recharge events, as supported by both the CCFs and graphical comparison between Q and EC (Fig. 11). However, the analyses of spring discharge variations and the VESPA index highlighted some differences between Nanto and Mossano springs. In fact, Mossano spring is composed of three systems with different degrees of permeability, whereas the Nanto spring reservoir is characterized by two systems. Another difference is that karst behavior prevails in Mossano spring (the conduit and intermediate systems), while fracture/matrix behavior prevails in Nanto spring. This result is further supported by the VESPA index, which delineates Nanto spring as a type C drainage system and Mossano as a type A drainage system.
The results obtained can be extended to the other perennial springs affecting the south-eastern part of the Berici Mountains in which the geological and hydrogeological conditions are the same as the Nanto and Mossano karst systems. However, a monitoring program is indispensable in order to correctly establish the hydrogeological features affecting these perennial springs.
Although, in this work, 4 years of data of continuously monitored parameters were collected, only the continuous data from 1 year were indispensable to perform all of the analyses included in the proposed workflow.
The combined use of such methodologies allows for both support of the detail of the hydrogeological conceptual model and limitation of the loss of information or erroneous interpretation that can result from using only one technique. For example, the results achieved by the ρ coefficient do not allow the identification of the existence of the karst conduits for Nanto spring, delineating a type-C behavior ascribable to a matrix-fracture system. However, although the RCA, FDC and ACF individuated the existence of both karst conduits and the matrix-fracture system, the VESPA index approach allows one to confirm that the matrix-fracture system prevails over the conduit system.