Guillaume Cinkus
email: guillaume.cinkus@umontpellier.fr

KARST: Service National d'Observation du Karst PNRGC

Keywords: karst, hydrological functioning, classification, typology, time series analysis, multivariate analyses. SNO

Classification is a first-line tool for understanding the main characteristics of a natural system's response. We propose a new classification of karst systems hydrological functioning that is based on spring discharge time series and takes profit of spring discharge databases to encompass the high diversity of karst hydrological functioning. It discriminates six different classes based on three relevant indicators of karst hydrological functioning.

A core dataset made of 10 karst systems was first considered for the set-up of the classification. The spring discharge time series were investigated according to recession curves, statistical and signal analyses to identify relevant indicators of hydrological functioning. The selection of the most relevant indicators and the proposal of the classification were based on multivariate analyses. The classification was then tested on spring discharge time series of 78 karst systems located worldwide.

All the systems homogeneously spread among the six proposed classes, which highlights the relevance of the approach and the representativeness of the various classes of hydrological functioning. Results from the proposed methodology were finally discussed to explore its limitations and define guidelines for its application.

Identification of relevant indicators for the assessment of karst systems hydrological functioning: proposal of a new classification

Guillaume Cinkus 1 , Naomi Mazzilli 2 and Hervé Jourde 1 1 HydroSciences Montpellier (HSM), Université de Montpellier, CNRS, IRD, 34090 Montpellier, France 2 UMR 1114 EMMAH (AU-INRAE), Université d'Avignon, 84000 Avignon, France

Introduction

10% of the world's population is dependent on karst water resources for drinking water [START_REF] Stevanović | Global distribution and use of water from karst aquifers[END_REF]. Karst systems are underground entities that drain recharge water over a catchment towards a main outlet. The water is drained through conduits, fractures and matrix, which originate from the dissolution of the calcite deposits by acidic water from the surface. Understanding the functioning of these complex and heterogeneous systems is therefore a major challenge for long-term water resource management. Over the past century, different methods have been developed to analyse hydrological time series, and subsequently characterize the functioning of karst systems. These methods can be considered as a preliminary step in the development and design of hydrological models of karst functioning for sustainable water resource management.

Classification has been widely used in surface hydrology to characterize hydrosystems. Although the threedimensional properties of aquifers generate an additional complexity, some authors proposed different methodologies to classify them based on geological, morphological and hydrological functioning analyses [START_REF] Heath | Classification of Ground-Water Systems of the United States[END_REF][START_REF] Dahl | Review of classification systems and new multi-scale typology of groundwater-surface water interaction[END_REF][START_REF] Heudorfer | Index-Based Characterization and Quantification of Groundwater Dynamics[END_REF]. In many cases however, these classifications fail to address the complexity of karst systems, which are strongly heterogeneous and correspond to a wide diversity of hydrological functioning. For these reasons, many authors worked on classifications specific to karst systems, either based on geological and morphological analyses [START_REF] Mylroie | Hydrologic classification of caves and karst[END_REF][START_REF] Waltham | Engineering classification of karst ground conditions[END_REF][START_REF] Jouves | Speleogenesis, geometry, and topology of caves: A quantitative study of 3D karst conduits[END_REF][START_REF] Veress | Karst Types and Their Karstification[END_REF], hydrological response analyses [START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF][START_REF] Mangin | Pour une meilleure connaissance des systèmes hydrologiques à partir des analyses corrélatoire et spectrale[END_REF][START_REF] Soulios | Contribution à l'étude des courbes de récession des sources karstiques: Exemples du pays Hellénique[END_REF][START_REF] Bonacci | Karst springs hydrographs as indicators of karst aquifers[END_REF][START_REF] Kullman | Nové metodické prístupy k riešeniu ochrany a ochranných pásiem zdrojov podzemných vôd v horninových prostrediach s krasovopuklinovou a puklinovou priepustnosťou[END_REF][START_REF] Flora | Hydrogeological characterization and discharge variability of springs in the middle Verde river watershed, Central Arizona[END_REF][START_REF] Springer | A comprehensive springs classification system: Integrating geomorphic, hydrogeochemical, and ecological critera[END_REF][START_REF] Rashed | Assessing degree of karstification: a new method of classifying karst aquifers[END_REF][START_REF] Malík | Use of recession-curve analysis for estimation of karstification degree and its application in assessing overflow/underflow conditions in closely spaced karstic springs[END_REF] or even karst groundwater microbiological analyses [START_REF] Sinreich | Microbiological monitoring and classification of karst springs[END_REF].

The spring discharge of a karst system is considered as a base information in karst hydrogeology. It results from the combination of flows from the different compartments of the hydrosystem (soil, epikarst, unsaturated and saturated zones). The widely increased development of karst spring discharge monitoring offers the opportunity to study the relations between discharge and hydrological functioning in depth. Some authors proposed classifications of karst systems based on the analysis of spring discharge time series, either with visual interpretations [START_REF] Soulios | Contribution à l'étude des courbes de récession des sources karstiques: Exemples du pays Hellénique[END_REF][START_REF] Bonacci | Karst springs hydrographs as indicators of karst aquifers[END_REF], by calculating indicators of functioning [START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF][START_REF] Mangin | Pour une meilleure connaissance des systèmes hydrologiques à partir des analyses corrélatoire et spectrale[END_REF][START_REF] Flora | Hydrogeological characterization and discharge variability of springs in the middle Verde river watershed, Central Arizona[END_REF][START_REF] Springer | A comprehensive springs classification system: Integrating geomorphic, hydrogeochemical, and ecological critera[END_REF][START_REF] Rashed | Assessing degree of karstification: a new method of classifying karst aquifers[END_REF] or by interpreting the parameters of recession models [START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF][START_REF] Kullman | Nové metodické prístupy k riešeniu ochrany a ochranných pásiem zdrojov podzemných vôd v horninových prostrediach s krasovopuklinovou a puklinovou priepustnosťou[END_REF][START_REF] Malík | Use of recession-curve analysis for estimation of karstification degree and its application in assessing overflow/underflow conditions in closely spaced karstic springs[END_REF].

However, the aforementioned classifications have been developed by analysing only few karst systems or without considering the high diversity of karst hydrological functioning. Therefore, diversity in karst systems physical properties and hydrological functioning is not fully considered, which impairs the relevance of these classifications and raises the need for a more generic approach. This work aims to provide a new classification of karst systems hydrological functioning with the following key features: (i) a clear methodological basis, (ii) the analysis of a wide diversity of karst systems representative of contrasted hydrodynamics behaviours, (iii) an approach being relevant worldwide and in a scarce-data context (i.e. sites where there is little knowledge of the system, or only discharge monitoring for a few years).

In this paper, we took advantage of the recent release of spring discharge time series in databases such as the French SNO KARST [START_REF] Jourde | SNO KARST: A French Network of Observatories for the Multidisciplinary Study of Critical Zone Processes in Karst Watersheds and Aquifers[END_REF] or the WoKaS database [START_REF] Olarinoye | Global karst springs hydrograph dataset for research and management of the world's fastest-flowing groundwater[END_REF] to propose a new classification of karst systems hydrological functioning. The typology describes a system where one single hydrodynamic response to precipitation impulse is expected. The aim of the classification is to characterize the hydrological functioning of a system, but not to decorrelate the factors that influence the functioning. The paper is organized as follows. In Section 2, we define the general characteristics of the karst systems considered in this study. Section 3 presents the various tools and analyses considered for the characterization of karst systems hydrological functioning. The most relevant indicators of karst hydrodynamics are identified and presented in Section 4, using the discharge time series of 10 well-known karst systems that cover a wide range of hydrological functioning. Section 5 is devoted to multivariate analyses that are considered for the proposal of the new classification described in Section 6. The discussion in Section 7 aims to evaluate the relevance of the proposed approach applied to 78 karst systems and to highlight some of its limitations. Section 8 gives the conclusions.

General characteristics of the karst systems considered in this study

This section presents the data we used to develop and test the classification, which involve in two different datasets: (i) a core dataset for the assessment of the most relevant indicators of functioning and the design of the classification, and (ii) a complementary dataset for assessing the most efficient recession model, testing the classification and identifying its strength and limitations.

Core dataset

To ensure the quality of the study and its relevance to the problem, we performed spring selection on the basis of three criteria: (i) quality of the hydrodynamic monitoring, which is function of time-step, instrumentation, measurement uncertainty and length of the time series, (ii) diversity of the hydrological functioning among the karst systems, meaning that the final dataset should cover a wide range of hydrological functioning (e.g. related to dimensions of the catchment, rainfall, degree of karstification, hydrological functioning, etc.) and (iii) existing knowledge from prior studies, to ensure that the classification is consistent with the actual knowledge on system functioning. We selected 10 karst systems (Table 1), located in France. We retrieved data from several organizations: The French Karst National Observatory Service (SNO KARST), the Parc Naturel Régional des Grands Causses (PNRGC), Suez, and the DREAL of Bourgogne Franche-Comté. Selected systems have been the subject of several comprehensive studies based of methods such as geology, cartography, field observations, tracing, geochemistry, time series analysis and modelling.

Table 1: General characteristics of the selected karst system [START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF][START_REF] Moussu | Prise en compte du fonctionnement hydrodynamique dans la modélisation pluie débit des systèmes karstiques[END_REF][START_REF] Blavoux | Bilan, fonctionnement et protection du système karstique de la Fontaine de Vaucluse (sudest de la France)[END_REF][START_REF] Bakalowicz | Etude hydrogéologique de la source de l'Espérelle et de son bassin d'alimentation[END_REF][START_REF] Maréchal | Fonctionnement hydrogéologique du système karstique de la Fontaine de Nîmes en crue[END_REF][START_REF] Cholet | Fonctionnement hydrogéologique et processus de transport dans les aquifères karstiques du massif du Jura[END_REF][START_REF] Lorette | Groundwater-flow characterization in a multilayered karst aquifer on the edge of a sedimentary basin in western France[END_REF] and their associated discharge time series.

Complementary dataset

In order to check the relevance of the method and its capacity to differentiate karst systems functioning, the classification resulting from the analysis of the aforementioned well-known karst systems was tested on a complementary dataset of 68 karst systems with different characteristics (e.g. dimensions of the catchment, meteorological regime, climate, and karstification degree). We worked with springs discharge time series of 23 French karst systems coming from a database provided by the French state (Banque Hydro) and took the other 45 from the WoKaS (World Karst Spring hydrograph) database, which provides details of over 400 karst systems worldwide [START_REF] Olarinoye | Global karst springs hydrograph dataset for research and management of the world's fastest-flowing groundwater[END_REF].

The quality of these 68 springs discharge time series is appreciated according to low (C2) to very good (A) quality discharge data (Figure 1C) as proposed in [START_REF] Olarinoye | Global karst springs hydrograph dataset for research and management of the world's fastest-flowing groundwater[END_REF]. The dataset is fairly well distributed throughout the world, with karst springs discharge from 17 countries (Figure 1B). The considered karst systems are located in various climatic conditions. According to Köppen-Geiger classification [START_REF] Peel | Updated world map of the Köppen-Geiger climate classification[END_REF], these climatic conditions correspond to 12 different climates, the temperate oceanic (Cfb) being the most represented (63.2%, Figure 1A).

Presentation of the selected discharge time series analysis methods for the characterization of karst systems functioning

The discharge is directly related to the recharge and the emptying of the capacitive function of a system, but also depend on the system geometry and hydrodynamic properties [START_REF] Malík | Evaluating Discharge Regimes of Karst Aquifer[END_REF]. This section provides details about four methods for analysing discharge time series and the indicators that can be retrieved from these methods to characterize the functioning of a karst system.

Statistical analyses

Statistical indicators of discharge time series provide basic information about the overall functioning of a system. The most common are the mean, minimum, maximum, standard deviation, and the various quantiles.

The mean interannual discharge depends on both the dimensions of the catchment and the mean recharge; it can therefore be used to assess the dimensions of a system. The observed minimum and maximum discharges make it possible to understand the flow amplitude. However, significant uncertainty related to the extrapolation of extreme discharges at springs is associated with this indicator.

The comparison between various karst systems is facilitated when using nondimensionalized indicators, such as the coefficient of variation (CV), which is the ratio between the standard deviation and the mean of the discharge time series. [START_REF] Netopil | Ke Klasifikaci pramenu podle variability vydatnasti[END_REF] proposed the calculation of a "characteristic discharge", which is the ratio between the 0.9 quantile (discharge value that is exceeded 10% of the time) and the 0.1 quantile. This "characteristic discharge" was referred to as spring variability coefficient (SVC) in [START_REF] Stevanović | Karst Aquifers-Characterization and Engineering[END_REF]. The SVC is less affected by extreme values. The specific discharge (QS) corresponds to the ratio between the mean interannual discharge and the catchment area. It allows assessing karst hydrodynamics with an implicit consideration of the dimensions of the catchment. However, the area information is not always available and also depends on the quality of hydrogeological studies over the catchment of interest.

Recession curves analysis

The hydrograph of a flood recession corresponds to the period when the discharge gradually decreases as water is not replenished [START_REF] Toebes | On recession curves. Recession Equations[END_REF]. It is possible to distinguish two regimes: (i) the influenced (quickflow) regime, which corresponds to the period when the discharge is influenced by the rapid infiltration of water into conduits in the unsaturated zone; and (ii) the non-influenced (baseflow) regime, which begins when rapid infiltration ends, and corresponds to the draining of the saturated zone and less transmissive compartments of the system [START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF]. In the literature, analysis of recession curves is mainly used to estimate groundwater reserves [START_REF] Drogue | Analyse statistique des hydrogrammes de décrues des sources karstiques[END_REF][START_REF] Forkasiewicz | Régime de tarissement de la Foux-de-la-vis (Gard) étude préliminaire[END_REF][START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF], determine the hydrodynamic parameters of the aquifer [START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF], and provide information on flows, drainage, and karstification degree [START_REF] Drogue | Analyse statistique des hydrogrammes de décrues des sources karstiques[END_REF][START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF][START_REF] Kullman | Nové metodické prístupy k riešeniu ochrany a ochranných pásiem zdrojov podzemných vôd v horninových prostrediach s krasovopuklinovou a puklinovou priepustnosťou[END_REF][START_REF] Malik | Assessment of regional karstification degree and groundwater sensitivity to pollution using hydrograph analysis in the Velka Fatra Mountains, Slovakia[END_REF]Kresic, 2007;[START_REF] Malík | Use of recession-curve analysis for estimation of karstification degree and its application in assessing overflow/underflow conditions in closely spaced karstic springs[END_REF]. The karstification degree is an indicator of the karst maturity of a system, which is directly related to groundwater recharge, storage capacity, spring discharge dynamics and system connectivity.

Numerous models for the analysis of recession curves of discharge time series have been proposed (Table 2). [START_REF] Boussinesq | Essai sur la théorie des eaux courantes[END_REF] and [START_REF] Maillet | Essais d'hydraulique souterraine et fluviale[END_REF] made the first proposals. [START_REF] Horton | The Role of infiltration in the hydrologic cycle[END_REF] and [START_REF] Barnes | The structure of discharge-recession curves[END_REF] then developed recession models based on Maillet's equation. [START_REF] Coutagne | Étude générale des variations de débit en fonction des facteurs qui les conditionnent[END_REF] proposed an equation that simulates the discharge from a reservoir. [START_REF] Padilla | Relative Importance of Baseflow and Quickflow from Hydrographs of Karst Spring[END_REF] highlighted that Coutagne's equation struggles to reproduce the recession curves of karst systems, and thus introduced a new parameter Qc. This later parameter corresponds either to the discharge from poorly transmissive zones of the aquifer, or to the discharge from aquitards outside the karst system; it may also have no physical meaning. [START_REF] Drogue | Analyse statistique des hydrogrammes de décrues des sources karstiques[END_REF] proposed to approach the whole recession curve with a hyperbolic function. [START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF] developed a twoequation recession model, representative of the influenced and non-influenced regimes. [START_REF] Kullman | Nové metodické prístupy k riešeniu ochrany a ochranných pásiem zdrojov podzemných vôd v horninových prostrediach s krasovopuklinovou a puklinovou priepustnosťou[END_REF] approached recession curves by fitting a model based on a linear equation [START_REF] Kullman | Rezim podzemnych vod s turbulentnym prudenım v puklinovo-prasovom horninovom prostredı[END_REF] for influenced regime, and Maillet's equation for non-influenced regime. Based on Kullman work, [START_REF] Malík | Use of recession-curve analysis for estimation of karstification degree and its application in assessing overflow/underflow conditions in closely spaced karstic springs[END_REF] proposed a classification of karst systems functioning according to the number of equations required to model the recession and the values of the parameters. Other authors [START_REF] Samani | Analysis of spring hydrographs for hydrogeological evaluation of karst aquifer system[END_REF][START_REF] Griffiths | Streamflow recession in basins with multiple water storages[END_REF][START_REF] Ladouche | Système karstique des Cent Fonts[END_REF] proposed models inspired by the aforementioned models. We did not consider physical-based models as we do not have any information on reservoir geometry. Further details about empirical, physical-based models and recession curve analysis can be found in [START_REF] Fiorillo | The Recession of Spring Hydrographs, Focused on Karst Aquifers[END_REF].

In this study, we wanted to assess the variability of the hydrological response of the karst systems, which is only possible when accounting for different recession periods. For this reason, we dismissed the use of a Master Recession Curve (MRC) which aims to overcome the problem of recession variation by combining several recession curves into one and only.

The Mangin's model has been widely used to characterize karst systems, mainly because the author proposed a classification based on two indicators derived from the parameters of the model equation. The parameter α is assumed to characterize the draining of the capacitive function of the karst system, which corresponds in most cases to flow from the saturated zone. In case of a low karstification degree, α can be affected by flows occurring in the unsaturated zone, which may induce a lag in the response in the non-influenced regime [START_REF] Mudarra | Relative importance of the saturated and the unsaturated zones in the hydrogeological functioning of karst aquifers: The case of Alta Cadena (Southern Spain)[END_REF]. The indicator k is strongly linked to the α recession coefficient of the first component of the Mangin model (Table 2), which is applied only on the non-influenced regime. The indicator k is thus focused on the slow depletion of the aquifer and assumed to characterize the ability of a system to store and return recharge water. It is calculated with the following equation:

𝑘 = 𝑉 𝐷𝑌𝑁 𝑉 𝑎𝑛
Table 2: Summary of the main models developed to analyse recession curves, corresponding equations and comment on their most appropriate usage. Qt (L 3 T -1 ) is the discharge at time t (T), Q0 (L 3 T -1 ) the discharge at t=0 (T), α the recession coefficient (T -1 ), n (-) a constant, Qc (L 3 T -1 ) a constant discharge, QR0 (L 3 T -1 ) the baseflow extrapolated at t=0, q0 (L 3 T -1 ) the influenced discharge corresponding to the difference between Q0 and QR0, η (T -1 ) a constant characterizing the speed of infiltration (η = 1/ti, with ti (T) being the duration of the influenced stage), ε (T -1 ) a constant characterizing the concavity of the influenced part of the recession curve and β (T -1 ) a recession coefficient for the turbulent flow. L and T are the dimensions for the base quantities of length and time, respectively.

Model Equation Comment

Boussinesq (1903) 𝑄 𝑡 = 𝑄 0 (1 + 𝛼𝑡) 2 -Non-influenced stage -Surface water [START_REF] Maillet | Essais d'hydraulique souterraine et fluviale[END_REF] 𝑄 𝑡 = 𝑄 0 ⅇ -𝛼𝑡 -Non-influenced stage -Surface water [START_REF] Horton | The Role of infiltration in the hydrologic cycle[END_REF] 𝑄 𝑡 = 𝑄 0 ⅇ -𝛼𝑡 𝑛 -More suitable to surface water [START_REF] Barnes | The structure of discharge-recession curves[END_REF] 𝑄 𝑡 = ∑ 𝑄 0 ⅇ -𝛼 𝑖 𝑡 𝑘 𝑖=1 -More suitable to surface water [START_REF] Coutagne | Étude générale des variations de débit en fonction des facteurs qui les conditionnent[END_REF]

𝑄 𝑡 = 𝑄 0 [1 + (𝑛 -1)𝛼𝑡] 𝑛 (1-𝑛)
-Suitable for karst systems [START_REF] Padilla | Relative Importance of Baseflow and Quickflow from Hydrographs of Karst Spring[END_REF] 𝑄 𝑡 = (𝑄 0 -𝑄 𝑐 )[1 + (𝑛 -1)𝛼𝑡] 𝑛 (1-𝑛) + 𝑄 𝑐 -Suitable for karst systems -Qc strengthens Coutagne model Drogue (1972) 𝑄 𝑡 = 𝑄 0 (1 + 𝛼𝑡) 𝑛 -Suitable for karst systems [START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF] 𝑄 𝑡 = 𝑄 𝑅0 ⅇ -𝛼𝑡 + 𝑞 0 1 -𝜂𝑡 1 + 𝜀𝑡 -Suitable for karst systems -Associated classification [START_REF] Kullman | Nové metodické prístupy k riešeniu ochrany a ochranných pásiem zdrojov podzemných vôd v horninových prostrediach s krasovopuklinovou a puklinovou priepustnosťou[END_REF] 𝑄

𝑡 = ∑ 𝑄 0𝑖 ⅇ -𝛼 𝑖 𝑡 𝑘 𝑖=1 + ∑ ( 1 2 + |1 -𝛽 𝑗 𝑡| 2(1 -𝛽 𝑗 𝑡) ) 𝑄 0𝑗 (1 -𝛽 𝑗 𝑡) 𝑘 𝑗=1
-Suitable for karst systems -Associated classification

With VDYN the dynamic volume and Van the interannual mean yearly volume of water discharged at the spring. The dynamic volume is calculated by integrating the exponential function:

𝑉 𝐷𝑌𝑁 = ∫ 𝑄 𝑖 ⅇ -𝛼𝑡 𝑑𝑡 = 𝑄 𝑖 𝛼 ∞ 0
With Qi the discharge at the time ti (ti being the time at which the flow is considered to be laminar and also the beginning of the non-influenced regime) and α the recession coefficient. In his work, [START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF] suggested to characterize the capacity of dynamic storage with the maximum calculated VDYN, as it tends towards a stable value for a large number of analysed recession curves. The indicator i is used to characterize the capacity of a system to filter and attenuate the precipitation signal. It corresponds to the discharge attributed to the influenced regime (second component of the Mangin model, Table 2) two days after the flood peak. This discharge is expressed in relative proportion to the influenced discharge q0 and is calculated with the following equation:

𝑖 = 1 -2𝜂 1 + 2𝜀
With η a constant characterizing the speed of infiltration (η = 1/ti, with ti being the duration of the influenced stage) and ε a constant characterizing the concavity of the recession curve during the influenced stage.

The classification initially proposed by [START_REF] Kullman | Nové metodické prístupy k riešeniu ochrany a ochranných pásiem zdrojov podzemných vôd v horninových prostrediach s krasovopuklinovou a puklinovou priepustnosťou[END_REF] and updated by [START_REF] Malík | Use of recession-curve analysis for estimation of karstification degree and its application in assessing overflow/underflow conditions in closely spaced karstic springs[END_REF] differentiate systems by their karstification degree. The methodology consists to reproduce a recession curve by fitting one to several equations (either exponential or linear) and calibrate the α and β parameters of each formula. The karstification degree is then deduced from a table based on the presence of different flow sub-regimes (i.e. the number and nature of the equations) and the value of the α and β parameters. It ranges from 0.5 to 10 and is associated with a description of assumed structure of the system and karst groundwater circulation.

Correlational and spectral analyses

Correlational and spectral analyses are time series analyses that are used to study the frequency content of a signal (referred to as "simple analysis") and relations between signals (referred to as "cross-analyses") [START_REF] Massei | Investigating transport properties and turbidity dynamics of a karst aquifer using correlation, spectral, and wavelet analyses[END_REF]. The simple analysis consists of calculating the autocorrelation function of a signal and the corresponding spectrum (obtained using a Fourier transformation, the calculations are detailed in appendix A). The principle is to compare the signal with itself over an increasing time interval or shift [START_REF] Jeannin | Analysis of karst hydrodynamic behaviour using global approaches: a review[END_REF]. The cross-analyses examine the transformation of the input signal into an output signal [START_REF] Padilla | Study of hydrographs of karstic aquifers by means of correlation and cross-spectral analysis[END_REF]. Signal analyses, mainly developed by [START_REF] Jenkins | Spectral Analysis and its Applications[END_REF], [START_REF] Hannan | Multiple time series[END_REF], [START_REF] Brillinger | The Identification of Point Process Systems[END_REF] and [START_REF] Box | Time series analysis. Forecasting and control. Holden-Day[END_REF], were first applied to karst hydrology by [START_REF] Mangin | Pour une meilleure connaissance des systèmes hydrologiques à partir des analyses corrélatoire et spectrale[END_REF].

According to [START_REF] Mangin | Pour une meilleure connaissance des systèmes hydrologiques à partir des analyses corrélatoire et spectrale[END_REF], a karst system can be characterized by its response time to a unitary impulse (precipitation) and its inertia, which depends on both the volume of groundwater reserves and karstification degree of the karst system.

Simple correlational and spectral analyses allow determining three indicators of karst hydrological functioning: (i) the memory effect (ME), which is the shift k for an autocorrelation coefficient rk of 0.2. It translates variation in discharge over time, and is directly related to the inertia of the karst system [START_REF] Marsaud | Structure et fonctionnement de la zone noyée des karsts à partir des résultats expérimentaux[END_REF]; (ii) the regulation time (RT), which is the inverse of the bandwidth, i.e. the maximum ordinate of the spectrum divided by 2 (value of the integral of the function between 0 and +∞). It provides information on the duration of the influence of a unitary impulse [START_REF] Larocque | Contribution of correlation and spectral analyses to the regional study of a large karst aquifer (Charente, France)[END_REF][START_REF] Kovács | Geometry and hydraulic parameters of karst aquifers: A hydrodynamic modeling approach[END_REF], on the volume of groundwater reserves [START_REF] Marsaud | Structure et fonctionnement de la zone noyée des karsts à partir des résultats expérimentaux[END_REF], and makes it possible to assess the overall organization of flows in the system (e.g., conduits, fractures, and cracks) [START_REF] Jeannin | Analysis of karst hydrodynamic behaviour using global approaches: a review[END_REF]; (iii) the cut-off frequency or spectral band breadth (SBB), which corresponds to the frequency f at which the value of the spectrum sf becomes negligible. Beyond this frequency, the spectrum is equal to zero, and can be assimilated to noise [START_REF] Jeannin | Analysis of karst hydrodynamic behaviour using global approaches: a review[END_REF]. The cut-off frequency provides information on the ability of the system to filter unitary pulses [START_REF] Marsaud | Structure et fonctionnement de la zone noyée des karsts à partir des résultats expérimentaux[END_REF]. The results provide a general idea of how a karst system works, with an emphasis on the inertia of the system and its capacity to attenuate the recharge signal.

Analysis of classified discharges

The analysis of classified discharges provides information on flow regimes within a system, based on discharge monitoring at the outlet of a karst system [START_REF] Marsaud | Structure et fonctionnement de la zone noyée des karsts à partir des résultats expérimentaux[END_REF]. Based on empirical observations, [START_REF] Mangin | Etude des débits classés d'exutoires karstiques portant sur un cycle hydrologique[END_REF] suggested that the distribution of discharges (or logarithm of discharges) from karst springs can be approximated by a half-normal Gaussian distribution (the calculations are detailed in appendix B). He concluded that the comparison between quantiles of measured discharges and quantiles given by the halfnormal distribution should follow a straight line. According to this theory, any discontinuities of the line (corresponding to an inflexion point) indicate inhomogeneity in the functioning of the system, below or above a certain range of discharge. Such changes can occur at low or high discharges, and may be positive or negative. The interpretation proposed by Mangin is based on an extremely strong hypothesis, which is that the statistical half-normal Gaussian distribution properly describes the distribution of discharges from karst springs under a "homogeneous functioning".

The method allows identifying particular events inherent in karst hydrology, e.g. overflow at outlet, leakage to another system, storage and emptying phenomena, time varying extent of the recharge catchment. It also allows assessing the quality of the gauging station [START_REF] Grasso | Etude critique des méthodes d'analyse de la réponse globale des systèmes karstiques. Application au site de Bure[END_REF][START_REF] Marsaud | Structure et fonctionnement de la zone noyée des karsts à partir des résultats expérimentaux[END_REF][START_REF] Dörfliger | Guide méthodologique, Les outils de l'hydrogéologie karstique. Avec la collaboration de Ph[END_REF].

Analysis and selection of indicators of karst dynamics

Analysis of discharge time series were performed for the 10 karst systems aforementioned (core dataset) based on methods detailed in Section 3, with the goal to select the most relevant indicators of karst systems hydrodynamics.

Statistical indicators

The results of the statistical analyses highlight the diversity of hydrological functioning of the core systems (Table 3). The mean discharge allows to distinguish systems with a low discharge (Aliou, Baget, Fontainede-Nîmes, Mouline and Toulon), with a medium discharge (Durzon, Esperelle, Lods, and Mouthe), and with a large discharge (Fontaine-de-Vaucluse). Although the mean discharge is highly correlated with the dimensions of the catchment, it is also dependent of the precipitation and the hydrological functioning of the system. CV and SVC are highly correlated with a correlation coefficient R = 0.925 (p-value = 0.00012). Both can be related to the inertia of the system and allow to differentiate reactive systems (Aliou, Baget, Esperelle, Fontaine-de-Nîmes, Lods and Mouthe) from inertial systems (Durzon, Fontaine-de-Vaucluse, Mouline and Toulon). There is no evident relation between QS and the characteristics of the system. As the dimensions of the catchment are bypassed, we suppose that QS may be related to the karstification degree or the specific recharge (volume of water that goes to the aquifer, by unit area). As QS requires the knowledge about the area of the recharge catchment, which is either unknown or with high uncertainty, it is not retained as a relevant indicator for the classification.

Recession indicators

The dynamics of a karst system can be either (i) at infra-day time scale, meaning that the system is reactive with fast variations in discharge of the order of an hour, or (ii) at daily time scale, meaning that the system has high inertia and changes in discharge can be assessed on a daily basis. A comparison of results issued from recession curves analysis performed on hourly and daily time series showed that, even for reactive karst systems, the daily time scale provides sufficient accuracy.

We selected a total of 93 recession curves (Figure 2) from the overall time series of each system of the core dataset, with no distinction between seasonal or event-scale recessions, on the basis of the following conditions: (i) the peak flood discharge must be significantly high regarding the overall dynamics of the system. We suggest at least one tenth of the maximum discharge of the discharge time series, or greater than the interannual mean discharge. However, we do not exclude the eventuality of site-specific thresholds; (ii) there should be little or no disruption during the recession (e.g. precipitation leading to untimely peaks). In cases where the disruption was of short duration, data could be removed and replaced with a blank; and (iii) the recession must be complete, meaning it should include both the influenced regime and the entire noninfluenced regime (with some tolerance for high-inertia systems). In the specific case when karst spring behaviour is influenced by a particular hydrological functioning (e.g. the activation of an overflow outlet) that appears on the recession curve as a bending point, our selection concerned only the last, unaffected part of the curve, including the end of the recession limb. This approach ensured that the models, which are not mean to fit curves with bending points other than the one between the influenced and non-influenced regime, were correctly calibrated. Information loss was relatively minor, as discharges that were excluded from the analysis only represented a tiny part of the overall flow. The bending point required for Mangin's model [START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF] was defined manually and corresponds to the time t when the non-influenced regime begins (when flow is considered as part of the non-influenced regime).

Out of the 9 recession models presented in section 3.1.2, 4 models suit the study's requirements (we identified them as relevant for karst hydrodynamics analysis but also easy to automate). They correspond to the models from [START_REF] Drogue | Analyse statistique des hydrogrammes de décrues des sources karstiques[END_REF], [START_REF] Coutagne | Étude générale des variations de débit en fonction des facteurs qui les conditionnent[END_REF], [START_REF] Padilla | Relative Importance of Baseflow and Quickflow from Hydrographs of Karst Spring[END_REF] and [START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF], and will be further referred as Hyperbolic, Coutagne, Padilla and Mangin models, respectively. We tested the 4 models by examining their performance in fitting all 390 observed recession curves of both core dataset and complementary dataset (Figure 3), and by performing a sensitivity analysis. We selected Mangin model for the analysis of recession curves, as it provides a consistent and very good fit for all the recession curves (lowest NSE of 0.91), with a limited equifinality. Moreover, this model has been widely used for years and its indicators are well known by the community. Results show that these indicators clearly differentiate among systems of the core dataset (Figure 4). The analysis of i values, which allow assessing the capacity of a system to filter and attenuate the precipitation signal, revealed a relationship between the value of this indicator, and the saturation state of the system (Figure 5), which corresponds to the volume of water stored in both the saturated and unsaturated zones. Variability in i can be a consequence of either (i) the volume of water already stored in the karst system when the recharge occurs (i.e. the influence of the saturation state on the system connectivity), but also (ii) the variability of the overall organization of flow and groundwater storage between matrix, fracture and conduits in the different compartments of the system. We therefore decided to take account of variability in i values (or lack thereof) as an additional information of karst system hydrological functioning, by including a new indicator based on the variance of i.

Signal indicators

The correlational and spectral analyses were performed with a sampling step of 1 day and a maximum offset m of 125 days, according to the suggestions of [START_REF] Mangin | Pour une meilleure connaissance des systèmes hydrologiques à partir des analyses corrélatoire et spectrale[END_REF] on short-term and long-term analyses. The short-term analysis can be performed on systems that have at least one continuous year of daily discharge recording. As the 10 core systems are all in the same climatic context according to the Köppen-Geiger climate classification, the starting date for the discharge time series was set as the beginning of the hydrological cycle (1 September) to consider the seasonality. The results highlight the diversity of hydrological functioning of the core systems (Table 4). Aliou and Mouthe are very reactive systems that rapidly transfer a response proportional to the intensity and duration of unitary pulses (precipitation), with almost no dampening of the recharge signal. Because they exhibit a fast response to recharge events but are less reactive than the former, Baget, Fontaine-de-Nîmes and Lods are referred as low-inertia systems. Durzon, Esperelle and Mouline correspond to the category of mediuminertia systems, which are able to filter a greater-or-lesser proportion of unitary pulses, and dampen the recharge signal. This category encompasses a wide variety of hydrological functioning, as the medium RT value may translate either (i) a medium inertia or (ii) a high variability of hydrological functioning with both high inertia and low inertia responses resulting in an average RT. Fontaine-de-Vaucluse and Toulon are considered as high-inertia systems, as they have a high filtration capacity and noticeably dampen the recharge signal.

Indicators issued from the analysis of classified discharge

Within the core dataset (Table 4), this analysis hints that there is no change in the hydrological functioning of 5 systems (Aliou, Baget, Durzon, Esperelle and Mouline), but that there might be hydraulic or flow properties changes beyond a certain discharge for 4 systems (Fontaine-de-Nîmes, Fontaine-de-Vaucluse, Lods and Toulon). For 4 other systems (Fontaine-de-Nîmes, Lods, Mouthe and Toulon), there might be an overflow outlet, discharge to another system or a temporary storage of water within the karst system when the discharge reaches a certain value. These interpretations are confirmed in the literature for Fontaine-de-Nîmes [START_REF] Maréchal | Fonctionnement hydrogéologique du système karstique de la Fontaine de Nîmes en crue[END_REF], Fontaine-de-Vaucluse [START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF], Lods, Mouthe [START_REF] Cholet | Fonctionnement hydrogéologique et processus de transport dans les aquifères karstiques du massif du Jura[END_REF] and Toulon [START_REF] Lorette | Groundwater-flow characterization in a multilayered karst aquifer on the edge of a sedimentary basin in western France[END_REF].

This analysis, whose interpretation requires prior knowledge of the system or field observations to be fully relevant, does not seem appropriate for the classification of karst systems. The method is based on a strong hypothesis that may not be suitable for all systems. Although some relevant indicators of functioning can be retrieved from this analysis (e.g. activation of an overflow outlet, flow to another system, …), the interpretation is overly influenced by (i) the quality of the discharge-water height relationship, and (ii) the subjective vision of the operator who is performing the analysis. Moreover, when there is poor or no prior knowledge of the functioning of a system, it is very difficult to identify the specific functioning behind a bending point on the curve (e.g. the difference between a bending point due to the activation of an overflow outlet, or due to uncertainties on ungauged discharges). Thus, we choose to not include the indicators of functioning retrieved from this analysis in the proposed methodology for the classification of karst system hydrological functioning.

Multivariate analyses

The aim of this section is to gain insights into the dataset and the relations between indicators. We applied two unsupervised techniques on a dataset consisting of 9 variables and 10 observations. The observations correspond to the 10 core karst systems, and the variables are relevant quantitative indicators of functioning resulting from the application of the different methods of discharge time series analysis (excluding the indicators resulting from the analysis of classified discharge). The selected indicators are kmax, imean, IR, αmean, ME, RT, SBB, CV and SVC. Terms "max" and "mean" correspond to, respectively, maximum and average values of the indicator over all the recession analysed.

Principal component analysis

Principle

Principal component analysis (PCA) is a multivariate method that aims to reduce the dimensions of an observation space by producing "principal components", which are linear combinations of initial variables of a dataset that retains the most possible variation. Principal components are uncorrelated with each other (i.e. orthogonal to the previous one) and are ordered according to the amount of variance explained by the combination [START_REF] Everitt | An introduction to applied multivariate analysis with R[END_REF]. The results of a PCA can be seen as a small number of new variables that contains most of the information of a large number of initial variables. The interpretation of a principal component is realized by looking at the correlation between initial variables and the component, i.e. assessing the contribution of each variable. This analysis provides information on trends in the dataset and allows identifying eventual complementarity between initial variables.

Results

The results of the PCA are presented as a biplot of the first factorial plane (Figure 6), which explains 87.9% of the total variance of the dataset.

The first principal component (PC1) is strongly correlated with all indicators except IR: (i) kmax informs about the capacity of dynamic storage, (ii) imean, ME, RT and SBB are related to the capacity to attenuate the precipitation signal, (iii) αmean characterizes the draining dynamic of the capacitive function and (iv) CV and SVC are about flow dynamics. Together, these indicators can be seen as characterizing the inertia of a system, which would correspond to the first principal component (PC1, 73.4% of variance explained). The second principal component (PC2, 14.5% of variance explained) is mainly correlated with IR and thus would reflect the variability of the hydrological response of a system. The third principal component (PC3, 5.8% of variance explained) is most correlated with αmean, which characterizes the draining dynamic of the capacitive function of a system.

It is worth mentioning that the high correlation between several input parameters (e.g. CV and SVC) may have a substantial effect on the PCA results, particularly on the variance explained by the principal components (i.e. PC2 and PC3 could have had a higher variance), by overemphasizing the contribution of redundant indicators.

Systems from the core dataset are mostly scattered along PC1 with rather inertial systems on the positive part (e.g. Fontaine-de-Vaucluse, Toulon) and reactive systems on the negative part (e.g. Aliou, Mouthe), and to a lesser extent along PC2 with "highly variable hydrological response" systems on the positive part (e.g. Esperelle, Durzon, Mouline) and "more regular, steady hydrological response" systems on the negative part (e.g. Fontaine-de-Vaucluse, Toulon).

Clustering

Principle

The purpose of data clustering is to identify clusters that contain observations or objects with similar characteristics [START_REF] Jain | Data clustering: a review[END_REF][START_REF] Govender | Application of k-means and hierarchical clustering techniques for analysis of air pollution: A review (1980-2019[END_REF]. Clustering analysis can be used to identify archetypes, and offer a better understanding of the structure within a dataset [START_REF] Halkidi | On Clustering Validation Techniques[END_REF]. This technique is considered unsupervised, because it is not based on predefined classes or examples that would give an idea of the structure of the dataset (Berry and Linoff, 1996).

We selected a Ward hierarchical clustering method for performing the analysis, which consists of a succession of binary fusions that minimize between-cluster variance until one cluster remains [START_REF] Murtagh | Ward's Hierarchical Agglomerative Clustering Method: Which Algorithms Implement Ward's Criterion?[END_REF]. The Ward distance is equal to [START_REF] Tufféry | Data mining and statistics for decision making[END_REF]:

𝐷(𝐴, 𝐵) = 𝑑(𝑎, 𝑏) 2 𝑛 𝐴 -1 + 𝑛 𝐵 -1
With D the Ward distance between two clusters A and B that have centers of gravity a and b and frequencies nA and nB. The analysis was realized with standardized data and Euclidean distance as measure of dissimilarity, which is calculated with the following equation:

𝑑 𝑒𝑢𝑐 (𝑥, 𝑦) = [∑(𝑥 𝑗 -𝑦 𝑗 ) 2 𝑑 𝑗=1 ] 1 2
With d the distance between two points x and y of a d-dimensional dataset, and xi and yi the values of the j th attribute of x and y, respectively [START_REF] Gan | Data clustering. Theory, algorithms, and applications[END_REF].

This method is suitable for clusters of different sizes and shapes and provides a graphical representation (dendrogram) which helps for understanding the clusters structure and how they are connected. The main advantages over a non-hierarchical method are (i) that it is not necessary to know the number of clusters prior to the analysis, and (ii) that the results do not depend on the choice of initial clusters [START_REF] Tufféry | Data mining and statistics for decision making[END_REF].

A way to assess whether a variable is relevant to characterize a cluster is to realize a value-test or v-test [START_REF] Lebart | Multivariate descriptive statistical analysis. Correspondence analysis and related techniques for large matrices[END_REF]. For a given quantitative variable, it involves in comparing the mean x ̅ k for this variable of a particular cluster k to the overall mean X̅ , with the formula:

𝑣 = 𝑥 𝑘 ̅̅̅ -𝑋 ̅ √ 𝑠 2 𝑛 𝑘 ( 𝑁 -𝑛 𝑘 𝑁 -1 )
With v the result of the test, s² the overall standard deviation, N the total number of observations and nk the number of observations in the cluster k. A value of |v| greater than 1.96 corresponds to a p-value less than 0.05, which rejects the following hypothesis: the mean of the particular cluster is equal to the overall mean [START_REF] Lebart | Multivariate descriptive statistical analysis. Correspondence analysis and related techniques for large matrices[END_REF]. In this case, the variable is relevant to describe the group of observations in the cluster.

Results

Results are presented as a dendrogram. Four relevant clusters (corresponding to two majors clusters) have been identified (Figure 7). There is a great similarity between clusters and PCA results (Figure 6).

A v-test was performed to assess the indicators that best characterize each cluster (Appendix C, Table C1).

The major clusters A and B are differentiated on the basis of the following indicators: kmax, ME, CV, RT, imean, SBB, SVC and αmean (in order of importance), which are related to the capacity of dynamic storage and inertia of a system. According to the sign of the v-test results, cluster A includes reactive systems, with low to medium capacity of dynamic storage and cluster B includes inertial systems, with higher capacity of dynamic storage. The systems from cluster 1 are characterized by a high αmean value, corresponding to a fast draining of the capacitive function. None of the indicators clearly characterize the systems in cluster 2. The systems in cluster 3 are characterized by a high IR value, corresponding to a high variability of the hydrological response. The systems in cluster 4 are characterized by very high RT, ME and imean and a very low SBB, corresponding to a high attenuation capacity of the precipitation signal.

Confrontation of the results with the actual knowledge of the functioning of the systems

PCA and clustering gave similar insight into the functioning of the 10 core karst systems. Two major clusters were identified (A and B, Figure 7): (i) cluster A is characterized by systems with a highly reactive functioning, which can be divided into two sub clusters (1 and 2) by looking the draining dynamic of the function; (ii) cluster B is characterized by systems with inertial functioning, and can be divided into two sub clusters (3 and 4) by looking at the variability of the hydrological functioning (Figure 6, Figure 7).

The hydrological response of Aliou and Baget corresponds to well karstified systems [START_REF] Mangin | Pour une meilleure connaissance des systèmes hydrologiques à partir des analyses corrélatoire et spectrale[END_REF], which is consistent with their position in the major cluster A. The karst conduits network induces floods of short duration, that can have a high amplitude [START_REF] Mangin | Contribution à l'étude hydrodynamique des aquifères karstiques[END_REF], with a response time estimated to 7h after a precipitation event for Aliou and 14h for Baget [START_REF] Sivelle | Dynamics of the Flow Exchanges between Matrix and Conduits in Karstified Watersheds at Multiple Temporal Scales[END_REF]. [START_REF] Sivelle | Dynamics of the Flow Exchanges between Matrix and Conduits in Karstified Watersheds at Multiple Temporal Scales[END_REF] also found lower transfer coefficients for Baget in their reservoir modelling study. These results on the two systems support the difference between both clusters (1 for Aliou and 2 for Baget), as Aliou is characterized by a faster draining dynamic.

Cross-correlational analyses between precipitation and discharge performed by [START_REF] Cholet | Fonctionnement hydrogéologique et processus de transport dans les aquifères karstiques du massif du Jura[END_REF] on Lods and Mouthe systems showed low response times (17h and 10h, respectively) and high peak values (0.3 and 0.41, respectively), which are characteristics of a very reactive functioning. On Mouthe system, the higher draining dynamic of the capacitive function (cluster 1) is consistent with both the lower response time and higher peak value of the cross-correlation function than the ones obtained on Lods system (cluster 2).

Fontaine-de-Nîmes is a reactive karst system with a moderate degree of karstification [START_REF] Maréchal | Fonctionnement hydrogéologique du système karstique de la Fontaine de Nîmes en crue[END_REF]. The hydrological response is fast due to a high infiltration rate and a fast water transfer in a well-developed conduits network [START_REF] Maréchal | Karst flash flooding in a Mediterranean karst, the example of Fontaine de Nîmes[END_REF], which corresponds to the characteristics of cluster A. [START_REF] Fleury | Karst flash-flood forecasting in the city of Nîmes (southern France)[END_REF] found that the draining of the saturated zone was slow, which is consistent with the position of the system in cluster 2.

Esperelle karst system is described as significantly fractured [START_REF] Moussu | Prise en compte du fonctionnement hydrodynamique dans la modélisation pluie débit des systèmes karstiques[END_REF], and characterized by both a high impulse response height and dampened recession [START_REF] Pinault | Inverse modeling of the hydrological and the hydrochemical behavior of hydrosystems: Characterization of Karst System Functioning[END_REF]. This description is consistent with the one of cluster 2, which consists of reactive systems with moderate draining of the capacitive function.

Durzon, Mouline, Fontaine-de-Vaucluse and Toulon systems are included in the cluster B. Mouline is considered as a complex karst system, with long response times despite presence of flush flow effect that induces a quick transfer during winter [START_REF] Pinault | Inverse modeling of the hydrological and the hydrochemical behavior of hydrosystems: Characterization of Karst System Functioning[END_REF]. Although there is no explanation for quick transfer in Durzon system, we can suppose that its functioning is similar to Mouline as it is located in the same area. The existence of both slow and fast dynamics on these systems depending on recharge event and/or geometrical structure highlights a high variability of hydrological functioning, which is characteristic of cluster 3. Fontaine-de-Vaucluse and Toulon are differentiated in another cluster (cluster 4) due to the low variability of their hydrological functioning. The high inertia and homogeneous response of the Fontaine-Vaucluse system can be explained by the thickness of the non-saturated zone (800 m in average) and the large area of its catchment, estimated to be about 1160 km2 [START_REF] Ollivier | Challenges and Limitations of Karst Aquifer Vulnerability Mapping Based on the PaPRIKa Method-Application to a Large European Karst Aquifer (Fontaine de Vaucluse, France)[END_REF]). The Toulon system, defined as a complex karst system by [START_REF] Lorette | Groundwater-flow characterization in a multilayered karst aquifer on the edge of a sedimentary basin in western France[END_REF], is a multi-layer system with the discharge of an unconfined, fractured and karstified aquifer that is permanently fed by another confined aquifer. This continuous alimentation can explain the mostly homogeneous response of the system.

Proposal of a classification of karst systems hydrological functioning

In this section, we first present the classification based on the results of multivariate analyses. The classification is then applied on both core and complementary datasets to assess the relevance of the approach, regarding (i) the coherence with the well-known hydrological functioning of the karst systems in the core dataset, and (ii) the distribution of karst systems among the different classes from a worldwide perspective.

Classification of karst systems according to various types of hydrological functioning

The exploration realized with PCA and clustering confirmed the expected relations between the functioning of karst systems and indicator values. In both methods, karst systems are first differentiated from two main aspects: their capacity of dynamic storage and their capacity to attenuate the precipitation signal. The former is expressed through kmax and the latter through imean, RT, ME and SBB (n.b. the absolute values of the Pearson correlation coefficients of each pair are greater than 0.85, suggesting strong correlations). We chose the indicator RT to characterize the capacity to attenuate the precipitation signal (further referred as global inertia, as it considers the overall organization of flows in the system, the shape and the dimensions of the catchment and the saturation state of the system). Indeed, this indicator is more relevant than (i) imean, which is biased by the number of available recession curves, and (ii) ME and SBB because their assessments are somehow questionable as they rely on an arbitrary threshold and a subjective evaluation, respectively. A second element of differentiation between karst systems is the draining dynamic of the capacitive function of a system with the αmean indicator. In this case, the mean of α values seems relevant due to the relatively low amplitude of the values for a given system (Figure 4C). A third element of differentiation between systems is the variability of the hydrological response by quantifying the variability of i with the IR indicator. The proposed classification thus relies on the following three main characteristics of the hydrological functioning: (i) the capacity of dynamic storage, (ii) the draining dynamic of the capacitive function, and (iii) the variability of the hydrological response. These characteristics are assessed using (i) kmax, (ii) αmean, and (iii) IR indicators, respectively. Based on these 3 distinct characteristics of hydrodynamic behaviour, six classes are proposed to discriminate the hydrological functioning of karst systems (Table 5, Figure 8): (i) C1:

Poor capacity of dynamic storage, with fast draining of the capacitive function and substantial variability of hydrological functioning; (ii) C2: Poor capacity of dynamic storage, with fast draining of the capacitive function and low variability of hydrological functioning; (iii) C3: Poor capacity of dynamic storage, with moderate draining of the capacitive function and substantial variability of hydrological functioning; (iv) C4: Poor capacity of dynamic storage, with moderate draining of the capacitive function and low variability of hydrological functioning; (v) C5: Noticeable capacity of dynamic storage, with slow draining of the capacitive function and substantial variability of hydrological functioning; and (vi) C6: Noticeable capacity of dynamic storage, with slow draining of the capacitive function and low variability of hydrological functioning. We chose to not include the indicator RT in the classification methodology, as it seems that the global inertia of a karst system is relatively independent of its main characteristics of functioning, especially for systems in C3, C4, C5, and C6 (Figure 9B).

The proposed classification is based on recession parameters derived from Mangin's model for recession simulation. This model gave good simulations results over the 78 catchments considered in this study. However, the good performance of this model should be checked before any use of the classification over other systems. To gain insights into the functioning of such systems, we recommend the operators to make use of correlational and spectral analyses, or analysis of classified discharges, which can already give relevant information about global inertia and functioning of a system.

Application of the proposed methodology to 78 karst systems

The classification was applied on the complementary dataset presented in 2.2 and composed of 68 karst systems, plus the 10 core systems for a total of 78 karst systems (Figure 9A).

The most represented class is C3, which represents systems with a very low to medium capacity of dynamic storage, a moderate to slow draining of the capacitive function and a medium to high variability of the hydrological response. About 74% of the systems are characterized with a very low to medium capacity of dynamic storage (C1, C2, C3 and C4), 26% have a fast draining of their capacitive function (C1 and C2) and 67% have a medium to high variability of the hydrological response (C1, C3 and C5).

RT was calculated for each system (except for 1 system in C3 and 4 systems in C2 with too much gaps in the time series), with means of 19.5, 13.9, 31.8, 41.4, 43.8 and 46.3 days for, respectively, C1, C2, C3, C4, C5 and C6 (Figure 9B). The small increase of RT throughout the classes is consistent given the functioning description of each class (C1 being the most reactive and C6 the most inertial). The smaller mean in C2 is related to the structure of the classification: a high IR in C3 and C5 means that the hydrological response can be more reactive than expected, whereas in C1 it means that the hydrological response can be more inertial than expected. It results in a higher RT mean in C1 (over C2) as it includes systems with potential inertial responses. Finally, we observed that RT is biased for systems with long dry periods (Saint-Pierre or Lez systems, the latter being under anthropic influence), thus it is suggested to not calculate this indicator for these systems.

The systems are spread across all classes and types of hydrological functioning. It means that, even applied on a wider dataset, there is a relative representativeness of all classes and types of hydrological functioning.

The spread of the 10 core systems between the six classes somehow confirms the respect of the second criteria for spring selection (diversity of the hydrological functioning among the karst systems, see section 2.1) and the relevance of the proposed classification.

Discussion

The aim of this section is to take a step back on the methodology by questioning some of its and assessing their impact on the relevance of the classification.

Influence of the length of the time series on the classification

We performed a sensitivity analysis in order to assess the reliability of the classification regarding discharge time series of short length. For 9 systems of the 10 core systems (Toulon was not considered as there is only 5 years of monitoring), we defined 7 scenarios that range from 1-year (Y1) to 7-years (Y7) length of the discharge time series. The classification methodology was then applied on each scenario in order to compare the results to those obtained with the full-length (FL) discharge time series (Table 6).

The accuracy of each indicator increases with the length of the time series (Figure 10). kmax is the most consistent with a steady decrease throughout the years in the to the kmax value obtained for the fulllength time series (FL indicator). αmean becomes more stable and probably more reliable from Y5 even though it does not show any significant increase in precision for Y6 and Y7. IR is highly uncertain for the shorter time series (<Y4) but stabilizes from Y5 and become relevant for Y6 and Y7. Aliou, Baget, Durzon and Mouline were correctly classified from Y2 (Table 6), Esperelle and Mouthe from Y5 and Y6, respectively, while Fontaine-de-Nîmes, Fontaine-de-Vaucluse and Lods were not consistently classified at all. The results emphasize the benefits of long-term monitoring (as there are more chances of observing multiple flood events, as well as different meteorological conditions), but also highlight several limitations:

• The classification is not reliable if only one recession curve is considered, as IR would be • When the indicator is close to the threshold, it can mislead the classification (e.g. Esperelle and Fontaine-de-Nîmes, for which the variability in αmean can be explained by either a highly variable or a complex hydrological functioning, respectively);

• As the maximum observed discharge Qmax can change over time, this can affect the recession curves selection, since only recessions with a peak flood discharge of at least one tenth of the maximum observed discharge are considered for the analysis (e.g. Fontaine-de-Nîmes with Qmax that went from 8.2 to 16.5 m 3 .s -1 from Y1 to Y3); • As the mean interannual discharge Qmean may vary over years, this can induce an evolution of kmax value and thus modify the class of the karst system when the value is close to the threshold (e.g. Fontaine-de-Vaucluse, for which kmax oscillate around the 0.4 threshold with changes in Qmean during the years, despite being issued from the same recession curve). Based on these results, we suggest working with at least 3-years length discharge time series for the classification. These 3 years should be taken as a guideline and may differ notably depending on the system's dynamics. Indeed, very reactive systems may require only a few years for a definitive classification thanks to their high hydrodynamic variability, while the minimum length of the discharge time series required to reach satisfying classification may increase for very inertial systems. However, it remains appropriate to work with shorter time series when there is no alternative. Our analysis shows that 7 out of 9 systems are correctly or almost correctly classified at Y2 if we include those that are close to the threshold (Table 6).

Evaluation of the distance between a system and other classes

The uncertainties related to either (i) the length of the discharge time series, or (ii) the indicators that are close to the threshold, can be addressed by estimating the distance to the other classes. The distance Ds-C of a system s to an adjacent class C is measured with the following equation:

𝐷 𝑠-𝐶 = |𝐼 𝑡ℎ -𝐼 𝑐𝑎𝑙𝑐 | 𝐼 𝑡ℎ
Where Icalc corresponds to the calculated value of a given indicator and Ith corresponds to its threshold value.

The indicator to consider in the calculation is the one that is critical for the class definition (i.e. corresponding to the junction in the flowchart). A distance D lower or equal to 0.1 (10%) can be considered as close to the threshold. As the distance increase, the system is not likely to be related with the involved class. In relation to the structure of the classification, there is no need to calculate the distance between classes 1-5 and classes 2-6, as it is highly uncertain that a system has both a kmax higher than 0.4 and a αmean greater or equal than 0.03 (hatched area, Figure 11A). The formula allows for calculating the distance between (i) classes 1-2, 3-4 and 5-6 with IR, (ii) classes 1-3 and 2-4 with αmean, and (iii) classes 3-5 and 4-6 with kmax.

The distance between diagonal classes can be calculated by applying the Pythagorean Theorem:

𝐷 𝑠-𝐶 = √𝐷 𝑠-𝐶𝑥 2 + 𝐷 𝑠-𝐶𝑦 2
Where x and y correspond to two different classes. The formula allows for calculating distance between (i) classes 1-4 and 2-3 with αmean and IR (Figure 11C), and (ii) classes 3-6 and 4-5 with kmax and IR (Figure 11B). These notions of distances can be more easily appreciated on the Figure 12, which shows the 3-dimensional distribution of the karst systems in the kmax, αmean and IR space. To complement the discussion in section 7.1 related to the "close to the threshold" systems, we calculated the distance to other classes for the 10 core karst systems (Table 7). The results highlight that (i) Baget is close to C3, (ii) Fontaine-de-Nîmes is close to C4 and (iii) Fontaine-de-Vaucluse is close to C4. Regarding the 78 systems, only 12 systems are close to a threshold with a distance to another class lower or equal than 0.1. This distance measure allows reflecting on the classification results by highlighting potential threshold issues. Thus, we can tell if a system is clearly into a said class, or is close to one or two classes.

Beyond the classification

Correlational, spectral and classified discharges analyses can be performed post-classification to really exploit the maximum of information that can be obtained by analysing the discharge time series, in order to get a deeper knowledge of the functioning of a system.

We propose to use RT as an additional indicator to complement the above-described classification and to gain further insights into the global inertia of a karst system. We defined four ranges: (i) lower or equal than 15 days; (ii) greater than 15 days and lower or equal than 30 days; (iii) greater than 30 days and lower or equal than 45 days; and (iv) greater than 45 days, which are referred as (i) low, (ii) medium, (iii) high, and (iv) very high global inertia, respectively.

In a similar manner, we propose to perform the visual interpretation of the curve of classified discharges after the classification to have additional information about the presence or absence of major specific functioning.

For example, HR_0020 system has a kmax of 0.07, an αmean of 0.021 day -1 , and an IR of 0.34, corresponding to C3 (Appendix D). HR_0020 is thus described as a system with a very low to low capacity of dynamic storage, a fast draining of the capacitive function and a medium to high variability of the hydrological response. The regulation time of 15.5 days indicates that the system has a medium capacity to attenuate the precipitation signal (or medium global inertia). The bending point on the curve of classified discharges at 1.8 m 3 .s -1 may be due to the activation of an overflow outlet, a discharge to another system, or a temporary storage of water. The one at 30 m 3 .s -1 may translate the same processes, but it is not excluded that it can be related to uncertainties on ungauged discharges.

Conclusion

Our objective was to propose a new classification of karst systems hydrological functioning based on the analysis of spring discharge time series that is representative of a wide diversity of karst systems. Several methods were considered to identify indicators of functioning. Multivariate analyses allowed the identification of relevant indicators that allow distinguishing the analysed karst systems. Three indicators were identified as the most relevant indicators and thus retained to characterize karst systems and propose a classification into six different types of hydrological functioning. The hydrological functioning of karst systems was distinguished according to their capacity of dynamic storage, the draining dynamic of their capacitive function and the variability of their hydrological functioning. The classification can be completed with two additional analyses to characterize the global inertia and highlight the presence or absence of major specific functioning.

The challenge of developing a relevant classification was addressed (i) by considering a core dataset of karst systems with a high diversity of hydrological functioning, (ii) by selecting the most relevant indicators of hydrological functioning and proposing a classification based on multivariate analyses, and (iii) by testing the relevance of the classification on spring discharge time series of 78 karst systems located worldwide.

As the methodology requires only spring discharge time series, which is the most common monitored data, and gives relevant results with only few years of monitoring, the classification can be used in scarce-data contexts. It can thus be seen as a modern tool for the classification of karst systems hydrological functioning, which provides researchers and stakeholders with a first insight into karst system functioning based on accessible and straightforward analyses. We emphasize that the proposed typology first aims to describe the hydrological functioning of a system where one single hydrodynamic response to precipitation impulse is expected, but remains useful even when two or more responses are observed.

One perspective of this work would be to provide a database that allows both comparing the hydrological functioning of several karst systems and then proposing a link between the developed classification and the design of lumped parameter models. Research perspectives include the study of the relation between classification and the relevant structures and parameters of models for rainfall-discharge simulation.

Appendices

A.

Calculation details for the correlational and spectral analyses

The autocorrelation function rk is calculated with the following equation:

𝑟 𝑘 = 𝐶 𝑘 𝐶 0
With the autocovariance function Ck:

𝐶 𝑘 = 1 𝑛 ∑(𝑥 𝑖 -𝑥̅ )(𝑥 𝑖+𝑘 -𝑥̅ ) 𝑛-𝑘 1
With k the shift (0, 1, 2, …, m), n the length of the series, m the maximum shift possible (generally m < n/3), xi the i th element of the series, xi+t the (i+t) th element of the series and x ̅ the mean of the series. The correlogram correspond to the plot of rk versus k.

The spectrum sf is calculated with the following equation:

𝑠 𝑓 = 2 [1 + 2 ∑ 𝐷 𝑘 𝑟 𝑘 𝑐𝑜𝑠 (2𝜋𝑓𝑘) 𝑚 𝑘=1 ]
With f the frequency (f = j/2m for daily timestep), rk the autocorrelation function and Dk a weighting function (Tukey-Hanning window) to ensure that the estimated sf is not biaised [START_REF] Mangin | Pour une meilleure connaissance des systèmes hydrologiques à partir des analyses corrélatoire et spectrale[END_REF]:

𝐷 𝑘 = 1 + 𝑐𝑜𝑠 ( 𝜋𝑘 𝑚 ) 2
The spectrum is represented on a plot of sf versus f.

B. Calculation details for the analysis of classified discharges

The procedure to obtain the curve of classified discharges involves in (i) the quantiles calculation of the empirical distribution function, (ii) the calculation of the corresponding variable from the reference distribution function, (iii) the graphical representation of the relation between the quantiles of the empirical and reference distribution functions, and (iv) the choice of the x-axis scale, either arithmetic or logarithmic (corresponding to a normal or log-normal adaptation of the reference distribution, respectively).

The repartition function corresponding to the cumulative probability density regarding the standard normal distribution is:

𝑃(𝑋 ≤ 𝑧) = 1 2 [1 + ⅇ𝑟𝑓 ( 𝑧 √2 )]
For a half-Gaussian distribution:

𝑃(𝑋 ≤ 𝑧) = ⅇ𝑟𝑓 ( 𝑧

√2 )

The observed discharges are plotted on the x-axis of the quantile-quantile plot.

C.

Results of the v-test applied on clusters A, B, 1, 2, 3 and 4.

Figure 1 :

 1 Figure 1 : Details about the complementary dataset, regarding (A) climate (Köppen-Geiger classification), (B) countries and (C) quality of the time series. (D) is the legend for the different abbreviations.

Figure 2 :

 2 Figure 2 : Selected recession curves for each of the 10 karst systems. n corresponds to the number of recession curves identified over the available discharge time series of each spring (a total of 93 recession curves is considered).

Figure 3 :

 3 Figure 3 : Comparison of the performance of the models with respect to observed discharge over all recessions of both core and complementary dataset: (A) Boxplot of relative errors for influenced and non-influenced flow regimes. (B) Boxplot of Nash-Sutcliffe Effiency (NSE).

Figure 4 :

 4 Figure 4 : Boxplot of the values of the indicators proposed by Mangin (1975) after analysing the spring discharge recession curves for each of the core systems: (A) i, (B) k, and (C) α.

Figure 5 :

 5 Figure 5 : Variability in i for the core systems. The standardized ratio between the mean discharge of the month preceding the flood and the mean interannual discharge (Q1M) is used as a proxy of the saturation state of the system.

Figure 6 :

 6 Figure 6 : Biplot of the first two principal components PC1 and PC2 resulting from the PCA performed on the dataset (9 variables for 10 observations). The points correspond to the PC1 and PC2 scores of each observation (i.e. karst systems). The arrows represent the correlation of the variables (i.e. indicators of functioning) with PC1 and PC2. The brown circle indicates the theoretical maximum extent of the arrows. Clusters (colours) are referring to the results of the hierarchical clustering performed on the dataset (see section 5.2).

Figure 7 :

 7 Figure 7 : Dendrogram resulting from the Ward's hierarchical clustering performed on the dataset (9 variables for 10 observations). The X-axis gives the distance between observations and/or clusters. The nodes (or vertical bars) indicate at which value two observations and/or clusters are merged. Four relevant clusters (1, 2, 3 and 4) are identified with different colours, and "major" clusters A and B correspond to the merging of clusters 1-2 and clusters 3-4, respectively.

Figure 8 :

 8 Figure 8 : Flowchart for the classification of karst systems based on three indicators issued from the recession curves analysis.

Figure 9 :

 9 Figure 9 : Results of the classification for the systems in the core and complementary datasets. (A) The 6 classes are depicted on the x-axis and the number of systems by class on the y-axis. The colour is related to the capacity of dynamic storage, the outline is related to the draining dynamic of the capacitive function and the pattern is related to the variability of the hydrological response. (B) On the right-upper side, the boxplot shows the distribution of RT among the systems.

Figure 10 :

 10 Figure 10 : Boxplot of the mean relative error of each indicator across 9 core systems and for the different scenarios, regarding the indicators issued from the analysis of the full-length discharge time series. Each indicator on the x-axis is associated with 7 boxes that are sorted from shorter (left) to wider (right) length.

Figure 11 :

 11 Figure 11 : Distribution of the 78 systems according to each pair of indicators: (A) αmean vs. kmax, (B) IR vs. kmax, and (C) IR vs. αmean. The colours correspond to the different classes.

Figure 12 :

 12 Figure 12 : 3D distribution of the 78 systems according to each indicator. The associated rectangles show the boundaries of each class. The colours correspond to the different classes.

Table 3 :

 3 Results of the statistical analyses for the core systems.

	System	Mean discharge	Minimum discharge	Maximum discharge	CV	SVC	QS
		(m 3 .s -1 )	(m 3 .s -1 )	(m 3 .s -1 )	(%)		(mm.d -1 )
	Aliou	0.45	0	28.91	190.4	32.1	3.20
	Baget	0.49	0.02	10.10	147.2	14.0	3.22
	Durzon	1.63	0.72	16.33	59.8	3.3	1.21
	Esperelle	1.11	0.16	14.50	139.8	11.8	1.05
	F-de-Nîmes (Fdn)	0.54	0	16.52	221.9	40.4	1.04
	F-de-Vaucluse (Fdv)	17.54	2.79	85	71.3	5.9	1.36
	Lods	1.01	0.15	8.54	127.2	19.6	2.48
	Mouline	0.51	0.19	3.15	45.7	2.6	1.37
	Mouthe	1.92	0.01	18.28	120.1	24.4	3.32
	Toulon	0.46	0.27	0.99	37.6	2.5	0.40

Table 4 :

 4 Results of correlational and spectral simple analyses and interpretation of the classified discharges curve for the core systems. (A) Systems with no apparent, specific functioning; (B) systems in which the hydraulic or flow properties change beyond a certain discharge; (C) systems in which there is an activation of an overflow outlet, a discharge to another system, or a temporary storage of water that occur above a certain discharge.

	System	(day) ME	(day) RT	(day -1 ) SBB	Interpretation of the classified discharges curve
	Aliou	4.6	11.2	0.41	A
	Baget	17.6	24.4	0.34	A
	Durzon	49.9	41.4	0.24	A
	Esperelle	28.3	30.0	0.25	A
	F-de-Nîmes	18.9	17.8	0.33	B
	F-de-Vaucluse	81.4	67.8	0.13	B
	Lods	13.0	23.5	0.36	B, C
	Mouline	57.8	44.6	0.265	A
	Mouthe	7.9	11.7	0.415	C
	Toulon	101.8	86.1	0.08	B, C

Table 5 :

 5 Characterization of the karst systems hydrological functioning for the six defined classes, in terms of capacity of dynamic storage, draining dynamic of the capacitive function, and variability of the hydrological response, with the corresponding indicator values.

	Class	Capacity of dynamic storage	Draining dynamic of the capacitive function	Variability of the hydrological response	kmax αmean	IR
	C1	Very low to medium	Fast	Medium to high	≤0.4 ≥0.03 ≥0.25
	C2	Very low to medium	Fast	Low to medium	≤0.4 ≥0.03 <0.25
	C3	Very low to medium	Moderate	Medium to high	≤0.4 <0.03 ≥0.25
	C4	Very low to medium	Moderate	Low to medium	≤0.4 <0.03 <0.25
	C5	Medium to high	Moderate to slow	Medium to high	>0.4 <0.03 ≥0.25
	C6	Medium to high	Moderate to slow	Low to medium	>0.4 <0.03 <0.25

Table 6 :

 6 Results of the classification performed on the different scenarios of discharge time series of different length (from 1 to 7 years, Y1 to Y7). For each scenario, there are details about the results of indicators kmax, αmean, IR, the corresponding class, and the number of recession curves considered in the analysis. The full-length (FL) column corresponds to the results on the classification on the whole discharge time series.

	System	Indicator	Y1	Y2	Y3	Y4	Y5	Y6	Y7	FL
		kmax	0.02	0.01	0.02	0.02	0.02	0.02	0.02	0.02
		αmean	0.127	0.127	0.115	0.090	0.084	0.083	0.086	0.067
	Aliou	IR	0.02	0.02	0.04	0.18	0.18	0.18	0.18	0.18
		Class	2	2	2	2	2	2	2	2
		R. num	2	2	3	6	7	8	10	12
		kmax	0.06	0.06	0.06	0.06	0.07	0.07	0.07	0.08
		αmean	0.032	0.028	0.026	0.026	0.028	0.028	0.027	0.021
	Baget	IR	/	0.07	0.07	0.07	0.09	0.09	0.13	0.23
		Class	/	4	4	4	4	4	4	4
		R. num	1	2	3	3	5	6	7	13
		kmax	0.54	0.63	0.68	0.72	0.75	0.76	0.75	0.78
		αmean	0.005	0.004	0.004	0.004	0.004	0.004	0.004	0.004
	Durzon	IR	/	0.51	0.68	0.68	0.68	0.68	0.68	0.70
		Class	/	5	5	5	5	5	5	5
		R. num	1	3	5	7	8	9	10	10
		kmax	0.14	0.09	0.10	0.11	0.12	0.12	0.12	0.13
		αmean	0.030	0.034	0.034	0.034	0.029	0.029	0.027	0.022
	Esperelle	IR	0.55	0.57	0.57	0.57	0.57	0.57	0.57	0.62
		Class	3	1	1	1	3	3	3	3
		R. num	3	5	7	7	9	9	10	11
		kmax	0.01	0.06	0.04	0.04	0.05	0.05	0.05	0.09
	F-de-Nîmes	αmean IR Class	0.035 / /	0.023 0.36 3	0.035 0.21 2	0.035 0.21 2	0.031 0.25 1	0.031 0.25 1	0.031 0.25 1	0.021 0.25 3
		R. num	1	2	2	2	3	3	3	6
		kmax	0.38	0.35	0.37	0.41	0.35	0.35	0.37	0.40
	F-de-Vaucluse	αmean IR Class	0.003 / /	0.004 0.10 4	0.005 0.10 4	0.005 0.10 6	0.005 0.13 4	0.005 0.13 4	0.005 0.13 4	0.005 0.14 6
		R. num	1	2	3	4	5	6	7	9
		kmax	0.03	0.08	0.09	0.08	0.08	0.08	0.07	0.09
		αmean	0.029	0.018	0.018	0.016	0.021	0.020	0.021	0.021
	Lods	IR	/	0.02	0.02	0.18	0.18	0.22	0.22	0.28
		Class	/	4	4	4	4	4	4	3
		R. num	1	3	3	6	7	10	10	11
		kmax	0.58	0.49	0.54	0.55	0.57	0.56	0.56	0.59
		αmean	0.003	0.005	0.005	0.005	0.005	0.005	0.005	0.005
	Mouline	IR	/	0.43	0.58	0.59	0.59	0.59	0.59	0.60
		Class	/	5	5	5	5	5	5	5
		R. num	1	3	6	7	8	8	8	9
		kmax	0.01	0.02	0.02	0.02	0.02	0.02	0.02	0.02
		αmean	0.060	0.060	0.064	0.064	0.064	0.061	0.061	0.059
	Mouthe	IR	0.08	0.11	0.11	0.11	0.11	0.25	0.25	0.33
		Class	2	2	2	2	2	1	1	1
		R. num	3	5	6	6	6	8	8	10

Table 7 :

 7 Distance to the other classes for the 10 core karst systems.

	System	Class	C1	C2	Distance to C3 C4	C5	C6
	Aliou	C2	0.27	0	1.25	1.22	/	/
	Baget	C4	0.32	0.31	0.08	0	0.81	0.81
	Durzon	C5	/	/	0.95	2.04	0	1.81
	Esperelle	C3	0.27	1.52	0	1.50	0.68	1.64
	F-de-Nîmes	C3	0.30	0.30	0	0.01	0.76	0.76
	F-de-Vaucluse	C6	/	/	0.46 0.004	0.46	0
	Lods	C3	0.31	0.33	0	0.13	0.78	0.80
	Mouline	C5	/	/	0.47	1.49	0	1.42
	Mouthe	C1	0	0.31	0.98	1.03	/	/
	Toulon	C6	/	/	1.48	1.43	0.39	0

  Table C.1: Results of the v-test applied on clusters A, B, 1, 2, 3 and 4, with values of each indicators of functioning. Bold entries highlight values for which the p-value is less than 0.05. The sign of the v-test value indicates if the mean of the cluster is lower or greater than the overall mean.

				v-test value
	Indicator of functioning
		Cluster A Cluster B Cluster 1 Cluster 2 Cluster 3 Cluster 4
	αmean	2.00	-2.00	2.65	-0.17	-1.18	-1.26
	kmax	-2.59	2.59	-1.28	-1.55	1.58	1.59
	imean	-2.33	2.33	-1.76	-0.89	0.63	2.22
	IR	-0.57	0.57	-0.67	-0.02	2.15	-1.45
	CV	2.47	-2.47	0.95	1.69	-1.53	-1.50
	SVC	2.23	-2.23	1.42	1.07	-1.44	-1.30
	ME	-2.54	2.54	-1.43	-1.37	0.71	2.41
	RT	-2.39	2.39	-1.48	-1.18	0.43	2.50
	SBB	2.26	-2.26	1.75	0.83	-0.40	-2.37

Acknowledgments

We thank the French Ministry of Higher Education and Research for the thesis scholarship of G. Cinkus as well as the European Commission for its support through the Partnership for Research and Innovation in the Mediterranean Area (PRIMA) programme under Horizon 2020 (KARMA project, grant agreement number 01DH19022A). For the data provided, the French Karst National Observatory Service (SNO KARST) and especially D. Labat and A. Probst for the data gathered by A. Mangin at the Baget and Aliou springs, the Parc Naturel Régional des Grands Causses (PNRGC), G. Lorette for data monitored by Suez, the DREAL Provence Alpes-Côtes d'Azur and the DREAL Bourgogne Franche-Comté are also acknowledged.

Abbreviations