Introduction

Karst aquifers are heterogenous media where most of the flow takes place through a conduits system that conveys groundwater to one or a limited number of springs. This transmissive component of the karst aquifer represents only a few percent of the karst formation. Although spatial organization is largely unknown, observations of karst drainage networks always show a hierarchically organized system that converges towards the outlet. Consequently, at the scale of the karst system (Mangin 1975), the monitoring of the karst spring allows one to obtain different proxies of mass and pressure transfers for characterization of the functioning of the karst system.

Many methods and statistical approaches have been proposed in the past to describe the karst spring hydrograph (Maillet 1905; Boussinesq 1904; Drogue 1972; Atkinson 1977; Schöeller 1967), and especially since the 1970s with the application and adaptation of time series analysis techniques (e.g. Jenkins and Watts 1968; Box and Jenkins 1976) to hydrology (e. g. Yevjevich 1972) or karst hydrology (Mangin 1984) purposes. There are a lot of recent studies that still acknowledge these works, even if other methods accounting for non-linearity or non-stationarity (Labat et al. 2000) have since been proposed. Besides, if the rainfall/discharge relationship is a transfer process of primary importance in karst hydrology, other relationships using water level measurements in wells (Pinault et al. 2005; Pinault and Schomburgk 2006; Charlier et al. 2015), caves or surface streams (Bailly-Comte et al. 2008; Salerno and Tartari 2009), but also continuous physico-chemical or water quality measurements (Massei et al. 2006; Valdes et al 2006; Bailly-Comte et al. 2011), can be used with time series analysis techniques to characterize specific pressure or mass transfers occurring within the karst system.

Today, some methods based on rainfall/discharge time series analyses and karst spring hydrograph description are commonly used in applied hydrogeology for the classification of karst aquifers, the quantification of the groundwater resource (Marsaud 1997; Bakalowicz 2005) or the application of vulnerability methods (Kavouri et al. 2011). All these methods were mainly developed and firstly applied by Mangin (1971, 1975, 1984), Padilla and Pulido-Bosch (1995), with some adaptation given by El-Hakim and Bakalowicz (2007). These studies used software provided in the late 1970s by the CNRS (French National Centre for Scientific Research) subterranean Laboratory of Moulis (France).

As a result, numerous graphical, statistical or signal processing methods have been developed for decades to interpret hydrological time series, but there is no simple and standardized tool that can be used for this purpose, which is necessary for a rigorous comparison of results between case studies. This Technical Note presents XLKarst, which has been developed to provide a simple and easy-to-use tool to process a selection of proven methods that characterize the functioning of karst systems. Some of these methods can be applied using a spreadsheet, a numerical code or a statistical analysis software package, but there has been a lack of easy-to-use tools that gather all these methods with the aim to standardize their application. The recent work of Cinkus et al. (2023) can however be mentioned in this regard with an application of some methods limited to the analysis of the discharge time series, while the developed XLKarst tool is also used to interpret the rainfall/discharge relationship (or any input/output relationship of a black-box system that can be considered as a linear and time-invariant system).

The need for standardized application of proven methods was highlighted in hydrology by Ladson et al. (2013) when reviewing the widely used Lyne and Hollick (1979) baseflow index. The same observation can be made in karst hydrogeology for the application of signal processing methods, the use of a probability plot to interpret breakpoints in the cumulative distribution function of karst spring flow, or to analyze recession and baseflow dynamics in karst spring hydrographs. In the framework of the GeoERA RESOURCE project, WP5 – CHAKA: Typology of karst aquifers and recommendations for their management (GeoERA 2023), 15 national geological survey institutes from Europe aimed to develop a joint and harmonized methodology for that purpose, and apply these methods on their own case studies. The XLKarst tool has been written in this context in an easy to use and widely shared environment (MS Excel). This article presents the different functionalities of this tool with an application to the Fontaine de Nîmes karst system, France, which is part of the KARST observatory network (SNO KARST 2019) initiative. The first two sections are devoted to the description of the methods that can be found in the “Time Series Analysis” and “Discharge” menus of the XLKarst tool, with a short description of the use of the results in karst hydrology. The last section shows an application to the Fontaine de Nîmes Karst system with daily data of rainfall and discharge over 22 years.

Time series analysis

Summary statistics

As a first step, XLKarst allows to compute summary statistics which are often used as broad descriptors of hydrologic time series. In addition to common statistics (mean, median, standard deviation), the following statistics and descriptors are provided:

  • The spring variability coefficient (SVC, Flora 2004), which is the ratio between percentile 0.1 and 0.9. It can be used to describe the responsiveness of karst springs and is less sensitive to spurious data than the ratio between the maximum and the minimum values;

  • The coefficient of variation (CV) has also been used to describe the responsiveness of karst springs (Mangin 1975; Flora 2004). This coefficient is computed as the ratio of the standard deviation of a time series to its mean value. It expresses the level of dispersion around the mean, whatever the frequencies.

  • The memory effect (Mangin 1984), in days, which corresponds to the lag for which the autocorrelation function falls below \(0.2\). The memory effect is automatically computed by XLKarst for lags lower than 125 days using the estimate of the autocorrelation function given by Jenkins and Watts (1968). Memory effects higher than 125 days are not computed, since their estimates are biased by the annual cycle of recharge dynamics.

  • The base flow index (BFI), which is computed as the ratio, in volume, of the baseflow to the total flow. A lot of baseflow separation methods exist in the literature to compute the BFI. Recently, Ladson et al. (2013) proposed a standard approach using the Lyne and Hollick (1979) digital filter. The approach used a reflection of the time series of 30 days to address “warm up” issues, and up to 9 passes of the digital filter with a filter coefficient ranging from 0.9 to 0.98. The result of the BFI is a function of the filter coefficient and the number of passes, which are respectively 0.91 and 3 as default in XLKarst. It is computed on the entire period covered by the time series, which should start and end at the beginning of a hydrological cycle. The daily baseflow time series is also reported on the plot showing the discharge time series;

  • The regulation time, \({T}_{\mathrm{reg}}\), in days, which is only computed for daily time series. \({T}_{\mathrm{reg}}\) corresponds to the half of the spectral density function for the \(0\) frequency, which is computed as the discrete Fourrier transform of the autocorrelation function using a Tuckey lag window (Jenkins and Watts 1968), with a truncation point \(m=125\) days of the correlation function. \({T}_{\mathrm{reg}}\) expresses the relative contribution of long term processes with periods higher than 2 × 125 = 250 days to the total variance of the time series.

  • The ratio \({\sigma }_{250}/\sigma (\mathrm{\%})\), which is a new parameter that can be used to assess \({T}_{\mathrm{reg}}\). XLKarst uses a moving average filter of 250 days to assess the standard deviation of long term processes (\({\sigma }_{250}\)). The result is divided by the standard deviation (\(\sigma\)) of the raw time series over the same period. It is expressed in % and can be related to \({T}_{\mathrm{reg}}\) following Eq. 1:

    $${T}_{\mathrm{reg}}\approx 125\times {\left(\frac{{\sigma }_{250}}{\sigma }\right)}^{2}$$
    (1)

    The ratio \({\sigma }_{250}/\sigma\) can be computed whatever the time step of the time series, assuming that this time step is short enough to capture most discharge fluctuations.

XLKarst also proposes to compute all these statistics for a set of discharge time series to get a table summary of these descriptors (referred as “multiple analyses”).

Correlation and spectral analyses

Correlation and spectral analyses are computed using the estimates given by Jenkins and Watts (1968). For univariate time series, XLKarst computes the autocorrelation function, the spectral density function and the power spectrum. For bivariate time series, XLKarst computes the 2 correlation functions and the crosscorrelation function, the 2 power spectra, the gain, the phase, the cross-spectrum and the squared coherency spectra. All the spectral analyses are based on the discrete Fourrier transform of the correlation functions using a Tuckey lag window (Jenkins and Watts 1968). The Tukey lag window is also known as the Tukey-Hanning window, and is equivalent to the Hanning window used in many data analysis software packages for its good spectral properties (e.g. Hearn and Metcalfe 1995). It is used to minimize sampling and truncation errors in the spectral estimates.

Jenkins and Watts (1968) provide a \(\left(1-\alpha \right)\left(\%\right)\) confidence limits that is used as a test for non-zero correlation on the plots of autocorrelation and crosscorrelation, with \(1-\alpha =95\%\) as default. Approximate \(\left(1-\alpha \right)(\%)\) confidence intervals are also computed for gain and phase estimates (Jenkins and Watts 1968).

Use of correlation and spectral analyses in karst hydrology

The crosscorrelation function is the data analysis tool for the identification of a transfer function (Box et al. 1994). Correlation and spectral analysis tools provided by XLKarst are thus neither restricted to the analysis of the rainfall/discharge relationships nor to karst system characterization. Other causal relationships can be described by these statistical tools at various time steps. For instance, water temperature datasets were used to assess groundwater transit time (Bailly-Comte et al. 2011), turbidity time series were used as output to understand the turbidity mechanisms in a karst system (Bouchaou et al. 2002; Amraoui et al. 2003), discharge time series were also used both as input and output to describe surface water groundwater interactions along a river reach between two gaging stations (Bailly-Comte et al. 2008) or between river flows and karst spring flows (Larocque et al. 1998). Diffusivity of karst aquifers has also been assessed from the gain function using water level responses to earth tides and barometric pressure (Larocque et al. 1998; Marsaud et al. 1993). It was also used to describe the spatial variability of groundwater transfer in a karstic aquifer based on the rainfall/water level relationships in wells (Delbart et al. 2016), or for the analysis of various microclimates of karstic cave systems using air temperature, 222Rn activity, CO2 concentration and relative humidity (Bourges et al. 2006).

Both the memory effect and the regulation time have been used by many authors to characterize karst system functioning (see for instance Larocque et al. 1998, Padilla and Pulido-Bosch 1995, Mathevet et al. 2004, Fiorillo and Doglioni 2010, Lorette et al. 2018). They have been recently used for karst aquifer classification purposes (Cinkus et al. 2021). The crosscorrelogram is used in these works as an image of the impulse response of the karst system to rainfall events, which allows one to easily compare various karst systems (Padilla and Pulido-Bosch 1995; Crochet and Marsaud 1997). It however requires that:

  • the input signal, i.e. the rainfall time series, does not show any autocorrelation. The use of daily rainfall time series is often recommended since shorter time steps would induce a higher autocorrelation in the rainfall time series,

  • the karst system may be conceptualized as a linear system, which can be discussed using the squared coherency spectrum from cross-spectral analyses

Mangin (1984) used the spectral density function instead of the power spectrum for the analysis of the rainfall-discharge relationship of karst systems. This allows one to better compare the spectra that have different variances, but this approach cannot be used to compute the cross spectrum. Some studies (Marsaud et al. 1993; Padilla and Pulido-Bosch 1995; Larocque et al. 1998) however used these unitary spectra based on the correlation to compute the cross-spectrum, which may end up with erroneous interpretation of the gain function. As a result, spectral density functions are only computed in XLKarst for univariate analyses.

Karst spring hydrograph analysis (“Discharge” menu)

Probability plot for cumulative distribution function of hydrologic time series

Mangin (1975) proposed a method based on the cumulative distribution function of discharge time series of karst springs to show systematic hydrodynamic behavior that is used for identifying the existence of specific flow regimes within the karst system. The method requires the use of several hydrological cycles to increase the robustness of the approach, especially for high flow values. It consists of the representation of the discharge cumulative distribution function on a half-normal probability plot. A constant time step is required, which must be chosen according to the dynamics of the processes involved. The use of the log of the discharge is also recommended by Mangin (1975). This method aims to highlight some systematic changes shown by breaking points and slope changes on the linearized cumulative distribution of the time series.

XLKarst proposes to use either the half normal or the normal-probability plots, with or without log transformation. Indeed, this method can also be used with water level time series for which normal-probability plots are often more appropriate. The user has to define the size of the frequency class, and the lowest value to be taken into account to focus on data that are not influenced by the recession (Mangin 1975). Figure 1 shows some of the main cases and the associated interpretation of the results. These interpretations must be taken with cautious and checked with results from others methods, including field characterization for overflow mechanisms. For low flows and low frequencies, the first slope change may be compared to the discharge that is often observed at the beginning of baseflow recession, denoted \({Q'}_{0}\) in the next section.

Fig. 1
figure 1

Interpretation scheme of probability plots for the cumulative distribution function (cdf) of discharge time series (modified after Mangin (1971), Crochet and Marsaud (1997) and Hakoun et al. 2020)). \({\mathbf{F}}_{\mathcal{Q}}\left(Q\right)=\mathbf{P}(\mathcal{Q}\le Q)\) represents the cdf of the random variable \(\mathcal{Q}\) of daily discharge evaluated for a discharge Q, and α1, α2 and α3 are the slopes of the linear trends used to interpret the cdf

Full size image

Recession curves analysis

3.2.1 Method

Mangin (1975) proposed a method based on karst spring flood recession analysis to separately characterize the flow dynamics occurring after the flood peak, i.e. (i) when flood recession is still influenced by recharge processes and (ii) during the baseflow recession. The conceptual model that was used by Mangin (1975) allowed one to link these flow dynamics to the functioning of two karst compartments, the infiltration zone and the phreatic zone respectively. While the relationship between the results of karst spring recession analysis and the volume of water stored in the phreatic zone may be questionable depending on the karst system, and particularly for karst systems with thick vadose zone, it is nevertheless a method that allows quantification of the groundwater resource, regardless of the origin of the water in the karst system. As pointed out by El-Hakim and Bakalowicz (2007), this method appears as the easiest to use and the most appropriate for karst spring recession curve analysis.

The method is usually applied with daily discharge time series. As a preliminary step, it requires the selection of flood recessions. This step is often poorly documented, while it can be a time consuming step for long time series. In addition, one has to consider flood recessions that are long enough to capture the whole baseflow dynamics, with little influence of possible secondary flood peaks to capture the quick flow component. Some subjectivity may arise in the results of the recession curve analysis from this step. As a result, with the objective of a standard approach, an automatic procedure has been proposed in XLKarst for the selection of flood recession curves.

Automatic selection of flood recession

A new approach for automatic selection of flood recessions is proposed by XLKarst using the 3 following parameters:

  • “Min duration” defines the minimum duration in days of a flood recession. The default value is 90 days. As shown by Abirifard et al. (2022), the dynamical volume will be more reliable for aquifers that show a long dry season, which also means that only long flood recessions must be considered to derive baseflow recession parameters.

  • “Min peak” defines the minimum value of the discharge that can be considered as a flood peak. The default value is the percentile 0.9 computed on the whole discharge time series.

  • The “lag” parameter is an optional parameter used when the flood dynamics is controlled by long term variations while short term variations induce many fluctuations of the flow during the recession. By default, a flood peak is defined at time t when \({Q}_{(t)}>{Q}_{\mathrm{Min peak}}\) and \({Q}_{(t-\mathrm{lag})}< {Q}_{(t)}> {Q}_{(t+\mathrm{lag})}\) with \(\mathrm{lag}=1\) time step. This lag parameter can be increased if the automatic selection of flood peaks failed to capture the actual flood peak.

The lowest value between two selected peaks is used to find the end of the recession period at \(t={t}_{\mathrm{end}}\).

This approach allows to extract flood recessions that can then be used to fit the mathematical model given by Mangin (1975).

List of parameters derived from recession cure analysis

Following Mangin’s method, the total flow \(Q(t)\) given by each flood recession is modeled as the sum of the infiltration function \(\psi (t)\) and the baseflow function \(\phi \left(t\right)\) using Maillet’s law:

$$\psi \left(t\right)={q}_{0}\frac{1-{~}^{t}\!\left/ \!{~}_{{t}_{\mathrm{i}}}\right.}{1+\varepsilon t},\; t\le {t}_{\mathrm{i}}$$
(2)
$$\phi\left(t\right)=Q_{R0}e^{-\alpha t}=Q_{0}^{'}e^{-\alpha(t-t_{\mathrm i})}$$
(3)
$$\mathrm{With }\;{Q}_{0}={q}_{0}+{Q}_{R0}$$
(4)

where:

  • \(t\) is the relative time since the flood peak;

  • \({Q}_{0}\) (m3/s) is the flood peak at the beginning of the recession, i.e. at \(t=0\);

  • \({Q}_{R0}\) (m3/s) is the simulated discharge using \(\phi \left(t\right)\) at \(t=0\). It is a fictive discharge that conceptually represents the discharge coming from storage at the beginning of the recession;

  • \({Q'}_{0}\) is the measured (and the simulated) discharge at \(t={t}_{\mathrm{i}}\);

  • \({q}_{0}\) (m3/s) is the initial infiltration rate, computed as the difference between Q0 and QR0;

  • \({t}_{\mathrm{i}}\) (d) represents the duration of the infiltration after which the baseflow can be modeled by Maillet’s law \(\phi \left(t\right)\);

  • \(\alpha\) (d−1) is the recession coefficient

  • \(\varepsilon\) is the non-dimensional coefficient of heterogeneity used in the infiltration function \(\psi (t)\).

At \(t={t}_{\mathrm{i}}\), recharge processes are supposed to have a negligible influence on karst spring flow dynamics. The spring discharge is thus controlled by the hydrodynamics within the phreatic zone, which means that the flow dynamics represents the emptying of the phreatic zone. Mangin (1975) called the dynamical volume \({V}_{\mathrm{d}}\) of the karst system the volume of water that can be theoretically drained after \(t={t}_{\mathrm{i}}\). Based on the analysis on spring recessions over several (> 3) hydrological cycles, the maximum value of \({V}_{\mathrm{d}}\) is used by Mangin (1975) to get an estimate of the dynamical volume of the karst system.

Some authors compute this volume from the beginning of the flood recession using the fictive extrapolation \({Q}_{R0}\) of the discharge for \(t\le {t}_{\mathrm{i}}\) (see for instance El-Hakim and Bakalowicz (2007), Paiva and Cunha 2020 or Olarinoye et al. 2022). The evolution of the baseflow is however unknown during that time and may greatly deviate from an exponential extrapolation (Bailly-Comte et al. 2010; Kovács 2021). As a result, with the objective of a standard approach, XLKarst computes the dynamical volume starting from \(t={t}_{\mathrm{i}}\), as recommended by Mangin (1975):

$$V_{\mathrm d}=\int_{t_{\mathrm i}}^{+\infty}\phi\left(t\right)=\frac{Q_{0}^{'}}\alpha$$
(5)

The Fig. 2 illustrates how these parameters can be derived from a recession curve.

Fig. 2
figure 2

Parameters derived from the recession curve analysis following the method of Mangin (1975)

Full size image

Automatic calibration procedure

As pointed out by Olarinoye et al. (2022), there is a need for automated, robust and objective methods for extracting karst spring recession components but also for fitting the parameters of any mathematical model that reproduce the whole recession dynamics. The model of karst spring recession proposed by Mangin (1975) requires one to fit 3 parameters: \({t}_{\mathrm{i}}\), \(\alpha\) and \(\varepsilon\).The automated procedure that is proposed by XLKarst for optimizing the parameters of Mangin’s model is to find the set of parameters that allow the maximum number of points to be observed around the theoretical model, regardless of the deviations between the model and the observations for the points that deviate from the theoretical model. These parameters are fitted using 2 characteristic times:

  • \({t}_{\upvarepsilon }\), which defines the number of points for \({t\le t}_{\upvarepsilon }\) to take into account for the optimization of the \(\psi\) function (Eq. 2)

  • \({t}_{\mathrm{end}}\), which defines the number of points for \({t}_{\mathrm{i}}<t<{t}_{\mathrm{end}}\) to take into account the optimization of the \(\phi\) function (Eq. 3)

A first loop is used over \({t}_{\mathrm{i}}\) to find the value that gives the highest number \(n\) of measurements which satisfy Eq. 6, where \(x\) is the threshold value set to \(5\%\) by default, \(X(t)=\mathrm{log}\left(Q(t)\right)\), \(Y(t)=\phi (t)\), and \({r}^{2}\) is the Pearson correlation coefficient between \(X\) and \(Y\).

$$\left|\frac{Y(t)-X(t)}{X(t)}\right|/{r}^{2}<x$$
(6)

This procedure allows one to find the best alignment between observation and simulation while discarding discharge fluctuations due to small recharge events during the recession. \(\alpha\) is then computed as either (default) the slope in a logarithmic scale between \(({t}_{\mathrm{i}}, {Q'}_{0})\) and the last value of the recession given by \(({t}_{\mathrm{end}},{Q}_{\mathrm{end}})\), or the slope of the regression line computed with \(X\) for \({t}_{\mathrm{i}}<t<{t}_{\mathrm{end}}\). Knowing \(\alpha\) and \({t}_{\mathrm{i}}\), \({Q}_{R0}\) and \({q}_{0}\) are computed from Eqs. 3 and 4 respectively.

A second loop is used over \({t}_{\varepsilon }\) to find its optimal value, with \(3\le {t}_{\varepsilon }<{t}_{\mathrm{i}}\). XLKarst computes for each value of \({t}_{\varepsilon }\):

  • \(\varepsilon\) as the slope of the regression line given by the linearization of Eq. 2:

    $$\frac{{q}_{0}}{Q\left(t\right)-\phi (t)}\times \left(1-{~}^{t}\!\left/ \!{~}_{{t}_{\mathrm{i}}}\right.\right)-1=\varepsilon t$$
    (7)

  • the number \(n\) of points which satisfy Eq. 6 with \(X=Q\), \(Y=\psi +\phi\) and \({r}^{2}\) the Pearson coefficient given by the linearization (Eq. 7).

Finally, the value of \({t}_{\varepsilon }\) that gives the highest \(n\) is selected, with the corresponding value of \(\varepsilon\).

The Nash criteria (Nash and Sutcliffe 1970) is also provided by XLKarst using the log of the measured and simulated data to discuss the ability of the model to reproduce the data. The log transformation is used to focus on the low values of the baseflow dynamics. It is however not used for the automatic calibration procedure due to its high sensitivity to small recharge events that may occur during the recession.

Other methods have recently been proposed by Olarinoye et al. (2022) for the optimization of Mangin’s recession model. These methods require setting a parameter for the duration of the period influenced by infiltration, or assuming that the baseflow recession begins when there are no significant variations (low CV) in the flow time series. Also worth mentioning are approaches based on the calculation of a master recession curve that requires setting a flow threshold to isolate the last component of the recession curve (see for instance Posavec et al. 2006). A sensitive analysis of various parameters controlling the assessment of the dynamical volume from recession analyses has been proposed by Abirifard et al. in 2022. They use a synthetic numerical model to describe the effect of recharge dynamics, aquifer and karst network geometries, and pumping and heterogeneous hydrodynamic properties of a synthetic karst aquifer. Among the results, this work clearly demonstrates that the delayed recharge may greatly influence the assessment of the dynamical volume, which claims for the use of Maillet’s law after a given time following the flood peak, i.e. when the exponential recession dynamic is effective. This delay is a consequence of the recharge processes, but also of the dynamics of recharge, so it varies from recession to recession. The use of a constant delay or a given discharge to isolate the baseflow recession after the flood peak can thus induce a bias in the calibration of the baseflow recession model. Setting such parameters is not required in the method proposed in the XLKarst tool. However, regardless of the automatic method chosen, a manual check of the results is still strongly recommended, which is why XLKarst cannot be used to perform the recession curves analysis without showing the graph of each model fit. The user can also modify the automatic model adjustment by changing the values of \({t}_{\upvarepsilon }\), \({t}_{\mathrm{i}}\), and \({t}_{\mathrm{end}}\). An example of automated model adjustment is shown later in the case study section.

Use of the recession curves analysis for karst system classification

Mangin (1975) uses the above recession curves analysis method to summarize information on recharge dynamics through the infiltration zone and on the flow regulation of the groundwater resource. Two parameters are used to classify karst systems accordingly:

  • The \(i\) parameter, called the “infiltration delay”, which is the mean value of \(\frac{1-{~}^{t}\!\left/ \!{~}_{{t}_{\mathrm{i}}}\right.}{1+\varepsilon t}\) for \(t=2\) days;

  • The \(k\) parameter, called the “regulating power”, which is the ratio of the maximum value of \({V}_{\mathrm{d}}\) to the mean annual volume drained by the karst spring. The assessment of \(k\) requires the use of long time series to ensure one captures the seasonal range in the assessment of the maximum value of \({V}_{\mathrm{d}}\)

The position of a given karst system in the \(i=f(k)\) diagram is used to infer some information on the karst development and to define the best methods for exploiting and managing the studied karst system (Mangin 1975; Bakalowicz 2005). Two classification diagrams are proposed by XLKarst: the original one from Mangin (1975) assuming \(k<0.5\) for karst systems, and the modified version given by El-Hakim and Bakalowicz (2007) where \(k\) is explicitly expressed in years in a log scale and can be higher than 1. XLKarst also provides uncertainties for these 2 parameters computed as follows:

  • An expanded uncertainty \({\delta }_{i}\) is given by Eq. 8, with \(n\) the number of \(i\) assessments and \(\mathrm{std}(i)\) the standard deviation of the \(i\) values:

$${\delta }_{i}=2\times \frac{\mathrm{std}(i)}{\sqrt{n}}$$
(8)

  • The \({\delta }_{k}\) uncertainty is computed as the error that would be done on the \({V}_{\mathrm{d}}\) assessment if the recession curve used to compute the value of \({V}_{\mathrm{d}}\) was discarded. A high uncertainty may reflect an inconsistent value for the dynamical volume that is used to compute \(k\). Such a result requires a closer inspection of the corresponding recession curve.

XLKarst also proposes a graphical thematic analysis to take into account an additional parameter through the size of the points. For instance, the use of the mean discharge of the karst spring can be useful to compare different karst systems in terms of groundwater resource availability. An example of karst system classification given by XLKarst is shown in the case study section.

Application to the Fontaine de Nîmes karst system

Site description and dataset

The Fontaine de Nîmes spring (FdN, identification No. BSS002ESPJ) is located in southern France in the center of the city of Nîmes (Fig. 3). It is the main outlet of a karst system developed in Cretaceous limestone of Hauterivian age. Its catchment area is about 55 km2 as defined by numerous tracing experiments (Fabre 1997) and a water budget calculation (Maréchal et al. 2005). The infiltration zone is about 40 to 80 m thick. Many speleological diving expeditions have been carried out from this spring since 1905 which makes it possible to know the organization of the karstic drainage in the downstream part of the system (Fabre 1997): the conduit network splits into 2 main branches, one toward the north-west and the other toward the north-east, resulting in western and eastern compartments of the catchment (Maréchal et al. 2008).

Fig. 3
figure 3

© IGN, 2023)

Fontaine de Nîmes (FdN) spring karst catchment. a Location map of FdN in the city of Nîmes, France. b Sketch of the full extent of the catchment. c Zoom on the monitoring network used in this study and the mapped karst network (Topographic map: Plan IGN v2

Full size image

Figure 3 shows the location of the well 9A (identification No. BSS002ESQL) that intersects the karstic drainage network. Unpublished speleological data show that the northeastern branch continues upstream towards this well. A tracer test was carried out in 2013 (Bailly-Comte et al. 2018) from this well intersecting a cave to better understand the hydraulic connection between the northeastern branch of the karstic network and the FdN spring. Water level measurements in well 9A combined with in-situ fluorescence measurements at the spring showed that the tracer was originally trapped in the karst drainage system at the injection point and was then mobilized after a major recharge event. Based on this tracer test result, an hydraulic connection between the well and the spring was observed for a discharge higher than 6 m3/s at the spring (Bailly-Comte et al. 2018). For higher flow conditions, Maréchal et al. (2008) show that the saturation of the fissured karst system induces many intermittent springs located in the whole karst basin.

The discharge time series stems from the conversion of water-level measurements at hourly or shorter time step in the Mazauric cave a few meters upstream from the Fontaine de Nîmes Sp. from September 1999 to September 2021. It is converted into discharge using a rating curve based on discharge measurements up to 16.3 m3/s with an in-situ ultrasonic flow meter at 15-min time steps from October 2004 to March 2005. A daily averaged discharge time series is used for the following application of XLKarst, with the daily rainfall provided by METEO-France at the NIMES-COURBESSAC station #30189001.

Results

Summary statistics

The summary statistics tool of XLKarst provides results for the Fontaine de Nîmes (FdN) discharge time series between 1st Sept 1999 and 31st Aug 2021 (as shown in Table 1 and Fig. 4).

Table 1 Statistical summary of the Fontaine de Nîmes (FdN) daily discharge time series from 1999 to 2021 (22 years)
Full size table
Fig. 4
figure 4

Plot of the daily discharge time series of the Fontaine de Nîmes and the associated baseflow component computed using the Lyne and Hollick (1979) filter (filter coefficient of 0.91 with three passes). Inset shows a zoom of the graph from September 2018 to March 2019

Full size image

Table 1 shows that the value of the SVC is higher than 20, which is very high and characterizes sharp flood peaks of short durations, relating to the flashy behavior of this karst system. Accordingly, the values of the memory effect and the regulation time are low (see the results of Padilla and Pulido-Bosch (1995) as reference). The σ250/σ parameter is 41%, this means that more than half of the standard deviation of the time series is due to short-term (period less than 250 days) fluctuations. The regulation time computed from Eq. 1 is 21 days, which is consistent with the value given by the spectral analysis (22 days, Table 1). The BFI is 0.42, which means that the volume of water that is stored in the karst system is relatively low, a significant part of the flow (58%) corresponding to the fast response to rainfall events. The daily computation of the Lyne and Hollick filter (Fig. 4) shows that only successive flood events produce a baseflow higher than 1 m3/s, which reaches 2 m3/s in December 2014 and December 2018 for flood peaks higher than 10 m3/s (Fig. 4), knowing that the maximum daily discharge is 16.51 m3/s. The combination of all these statistics and descriptors highlight the flashiness of the Fontaine de Nîmes karst system.

Correlation and spectral analyses

Figure 5 shows the results provided by XLKarst for the bivariate time series analysis using \(m=125\) days, corresponding to a short term analysis of the rainfall-discharge relationship of karst systems (Mangin 1984).

Fig. 5
figure 5

Results of the correlation and spectral analysis for the daily rainfall/discharge relationship for the FdN karst system. a Autocorrelation functions; b Crosscorrelation function; c and d Power spectra of the rainfall and discharge time series respectively; e Gain function; f Phase function; g Cross-spectrum; h Squared coherency function. The red dashed lines represent the 95% confidence intervals

Full size image

The black curve on Fig. 5a shows that the rainfall time series is not autocorrelated at a daily time step since the autocorrelation function \({r}_{\mathrm{xx}}\) falls to zero for small lags. The test of randomness (red dotted line) shows that the autocorrelation function may be considered as 0 for most lags > 1, except a peak above the confidence threshold occurring between 50 and 60 days. This information will be used in the following for the interpretation of the crosscorrelation (Fig. 5b). The corresponding power spectrum (Fig. 5c) shows that the variance of this time series is distributed all over the frequencies, even if it is not a perfectly flat spectrum that would characterize a pure random process. The daily rainfall can thus be interpreted as a random process for the analysis of the rainfall-discharge relationship, which is a key element for the interpretation of the crosscorrelation function. With a shorter time step, there would be a greater chance of observing a rainfall event if it was already raining at the previous time step.

In contrast, the daily discharge time series is clearly autocorrelated, with a memory effect of 15 days (grey curve on Fig. 5a). The test of randomness shows that the correlation is still significant for lags higher than 70 days. The corresponding power spectrum (Fig. 5d) shows that the variance of the daily discharge time series is mostly explained by a seasonal trend for a period higher than 250 days, i.e. a frequency lower than 4.0 10–3 cycle per day (cpd). It can also be used to compute a regulation time of 22 days, which supports that the daily time step is large enough to describe the main characteristics of the impulse response of this karst system.

The crosscorrelation function \({r}_{\mathrm{xy}}\) is clearly an odd function, which highlights the causal relationship between the rainfall and the discharge. The maximum of this function is 0.49 (Fig. 5b). It is reached for a lag of 1 day. For higher lags, it decreases rapidly to reach a non-significant value for a lag of 38 days. Compared to others karst systems, this crosscorrelogram characterizes a system with low inertia and a short response time, typical of a well-karstified system. This interpretation is consistent with the shape of the impulse responses that was used by Maréchal et al. (2008) to simulate the spring discharge using a lumped transfer function model. This means that this karst system is not able to regulate recharge events, which explains its high responsiveness to intense Mediterranean rainfall events, inducing karst flash-floods (Maréchal et al. 2008). Secondary peaks can be seen on the crosscorrelogram, especially for a lag around 55 days. These peaks are simply due to the non-random properties of the rainfall time series for these corresponding lags (black curve on Fig. 5a).

Gain (Fig. 5e) and cross-spectra (Fig. 5g) show a peak for the lowest frequencies. The infiltration zone, and in particular the evapotranspiration processes taking place there, act as a filter that amplifies seasonal effects. This effect appears on the spectrum for a period greater than 250 days, which corresponds to the long-term component of the spectrum at zero frequency. Conversely, short-term processes are filtered out by the infiltration zone, producing little or no response at source. The FdN karst system can thus be compared to a low-pass filter with a cut-off frequency close to 0.2 cpd.

The squared coherency spectrum (Fig. 5h) shows values that are lower than 0.6, even for low frequencies that carry most of the covariance. This means that the use of a simple transfer function that would transform the daily rainfall time series into daily discharge using a convolution integral would give poor results. In other words the daily rainfall-discharge relationship cannot be considered as a linear process, even if this assumption allows one to describe the main characteristics of the system through a correlative and spectral analysis.

Finally, the phase spectrum (Fig. 5f) mainly shows a linear trend according to frequencies, which is typical of a simple delay system. The group delay given from the linear slope of the phase function according to frequency is of 0.8 days with \({r}^{2}=0.97\). This result is consistent with the lag of 1 day given by the sharp peak of the crosscorrelation function.

This correlation and spectral analysis allows one to describe the impulse response, or transfer function of the system, with a mode at one day (which explains most of the transfer) and a length of 20 to 30 days (which controls the low inertia of this karst system).

Probability plot

Figure 6a shows the cumulative distribution function of the log of the daily discharge time series using a class interval of 0.5 m3/s and a minimum discharge of 0.5 m3/s. This latter value is relatively high to focus on hydraulic functioning that is not influenced by the baseflow recession. This distribution shows three linear segments. The two first ones characterize the hydraulic functioning of the FdN karst system for daily discharge lower than 13.5 m3/s, where a slight change on the distribution occurs for Q = 6 m3/s. Figure 6b shows an extract of water-level measurements in the northeastern branch of the karstic network compared to the FdN spring discharge at a hourly time step. This illustrates that the northeast branch of the karstic network is only active during flood conditions, and more precisely when the discharge at the FdN Sp. exceeds 6 m3/s. This induces a change in the drainage capacity of the FdN Sp, with lower head losses and a modification of the hydraulic head/discharge relationship. This effect induces a slight change on the slope of the hydrograph, and thus on the probability plot. This interpretation is also consistent with the tracer test result mentioned in the site description of Bailly-Comte et al. (2018).

Fig. 6
figure 6

a Cumulative distribution function of the daily discharge time series of the FdN karst spring plot on a half normal probability plot in log coordinates for class intervals of 0.5 m3/s and Q > 0.5 m3/s, the three red dashed lines show the two slope breaks on the distribution at 6 and 13.5 m3/s b Comparison of water level measurements in Well 9A and discharge measurements in FdN showing the activation of the northeast branch for a discharge higher than 6 m3/s

Full size image

A last break point is shown for discharge higher than 13.5 m3/s (red line no. 3 in Fig. 6a), which, according to the classification scheme (Fig. 1) is due to the whole saturation of the karst system and the activation of several overflow springs, as observed by Maréchal et al. (2008).

Recession curve analysis and classification of karst systems

The automatic selection of flood recessions using a “Min duration” of 60 days and a “Min peak” of 4 m3/s allows selecting 24 events (Fig. 7) for 22 hydrological cycles. A visual inspection in semi-log coordinates gives a global view on the whole dataset before editing each recession parameter. It also helps to identify some inconsistencies in the data for low flows. XLKarst can be used for this purpose by easily zooming in on specific events to check the consistency of the data before editing the recession parameters, and possibly removing a given event from the analysis.

Fig. 7
figure 7

Overall view of the discharge time series in semi-logarithmic axes showing the recession curve fits and their numbering

Full size image

Figure 8 shows an example of the results given by XLKarst with the automatic fitting procedure for the 16th recession curve on semi-log and linear plots.

Fig. 8
figure 8

a Semi-log and b linear plots of the recession curves showing the values of tε, ti and tend that are used to fit the recession model shown by the red line (automated procedure)

Full size image

The value of \({t}_{\varepsilon }\) is relatively short compared to \({t}_{\mathrm{i}}\) since the discharge measured between \({t}_{\varepsilon }=10\) days and \({t}_{\mathrm{i}}=34\) days shows some spikes due to a small recharge events. The fitting process allows one to reproduce the beginning of the recession, which is the critical period for the assessment of the \(i\) parameter. This recession analysis gives \(i=0.39\) with \(\varepsilon =0.71\) and \({t}_{\mathrm{i}}=34\) days.

The results given by the fitting of the 24 recession curves are shown by red dotted lines on Fig. 7 and can be easily extracted from XLKarst as a table that gives all the parameters of the recession curves analysis. All these values are used by XLKarst to compute the \(i\) and \(k\) parameters and their respective uncertainties: \(i=0.43 \pm 0.04\) and \(k=0.18 \pm 0.06\). These values are reported in the two classification systems proposed by XLKarst (Fig. 9) to compare the results with those given by Mangin (1975) to build up this classification scheme based on 5 karst systems. The size of the point is proportional to the log of the mean discharge of each karst system that ranges from 0.45 m3/s (Aliou) to 17.5 m3/s (Fontaine de Vaucluse) according to Mangin (1975).

Fig. 9
figure 9

Classification of karst systems based on the recession curves analysis a Classification of Mangin (1975), b the modified classification proposed by El-Hakim and Bakalowicz (2007)

Full size image

Some key interpretations of these diagrams are given by Mangin (1975), Marsaud (1997) and El-Hakim and Bakalowicz (2007). The position of the FdN karst system highlights its low regulating power which is consistent with statistical and time series analyses. However, the value of \(i\) is relatively high, which reflects some delayed infiltration through the infiltration zone and/or complex surface/groundwater interaction through ephemeral streams draining less permeable compartments. This value is significantly higher than the value of 0.18 proposed by Maréchal et al. 2005 based on 4 recession curves from 1998 to 2004. This new assessment is however more reliable considering the length of the time series that is used in this study.

Conclusion

XLKarst is a tool that provides a standardized approach for a simplified but rigorous application of some signal processing tools and recession curves analysis, which is necessary for a rigorous comparison of results between different case studies. Indeed, various interpretations of these methods can explain some discrepancies between studies which argues for the development of a standard procedure such as the XLKarst tool.

XLKarst addresses the main issues related to the characterization of the hydrodynamic behavior of a karst aquifer by gathering and standardizing various methods in a single and easy to use tool. For the spectral analysis, this tool clarifies the computation of gain and cross-spectra that were computed differently in previous studies. XLKarst also proposes a method for the selection of the recession curves and the calibration of the parameters, which tends to limit as much as possible the subjectivity of the results when applying the method of recession curves analysis proposed by Mangin (1975). Finally, XLKarst provides a standard method for calculating the dynamic volume that may have been also defined differently in other studies, with an estimate of its uncertainty.

This work also shows an application to the FdN karst system to highlight the use of the various methods provided by XLKarst with some key elements of interpretation. The 22 years daily time series of rainfall and discharge reveal that the FdN karst system is a poorly regulated system prone to flash flood events in a Mediterranean climatic context. This system has a complex hydraulic behavior, with a hydraulic connection to a temporary active karst drainage system and overflow processes during large floods.

XLKarst is an Excel VBA spreadsheet that can be freely downloaded on the XLKarst webpage hosted by the French Geological survey (BRGM 2021) at https://www.brgm.fr/en/software/xlkarst-excel-application-hydrodynamic-characterisation-karst-systems

This tool will evolve and be enriched with other features. The VBA code can be shared with people interested in this development upon request to the authors.