Quantitative classification of carbonate aquifers based on hydrodynamic behaviour

A quantitative classification of carbonate aquifers based on hydrodynamic behaviour is introduced. This type of classification is necessary to understand the physical functioning of carbonate hydrogeological systems and to provide a realistic interpretation of field data. Carbonate aquifers are generally considered as karst systems; however, geomorphology and aquifer geology alone are insufficient for determining hydrodynamic behaviour. Analysis of spring and well hydrographs based on analytical solutions is applied to establishing a quantitative classification. A base-flow recession coefficient is used as an indicator of hydrodynamic behaviour. Detailed numerical analyses suggest that carbonate systems can be classified into two distinct groups based on hydrodynamic behaviour. The physical processes depend on a combination of hydraulic and geometric parameters, and their functional relationships can be quantitatively determined. The proposed classification methodology involves making an assumption about aquifer type, estimating aquifer properties from hydrograph data, and comparing the results with field observations. The proposed classification methodology was applied to aquifers representing the two groups of carbonate systems. In both cases, the applied methods revealed crucial information about hydrodynamic functioning of the investigated systems. While the studied limestone aquifer showed karstic hydrodynamic behaviour, the investigation of a dolomite aquifer disproves a priori assumptions on karstic flow conditions. Dolomite aquifers represent an ambiguous group of carbonates and require caution in the selection of investigation tools and interpretation of hydrogeological data. The introduced methodology provides a reliable means of determining the hydrodynamic functioning of an aquifer and supports the quantitative classification of carbonate hydrogeological systems.


Introduction
Carbonate aquifers represent a special group of hydrogeological systems. Certain limestone rocks manifest intense karstification and develop characteristic karst features (sinkholes, shafts, cave systems, etc.) and typical karstic hydrogeological behaviour such as high flow velocities in karst conduits and abrupt changes in spring discharge. Other carbonate rocks seem to lack these features, and do not appear to manifest the erratic hydrogeological behaviour considered to be characteristic of karst systems (Mangin 1975;Zeizel et al. 1962;Klimchouk 2015). While the presence of karst geomorphological features usually indicates the existence of karstic groundwater flow (Kiraly 2003), the absence or invisibility of these features do not necessarily mean the lack of karstic groundwater flow conditions (Klimchouk 2015;Granger et al. 2001;Worthington 1999;Worthington et al. 2000).
Groundwater flow is influenced by the physical setting of the flow media. Geology influences hydraulic conductivity and porosity through the abundance, aperture size, orientation and connectivity of voids. These characteristics are determined by sedimentary processes, diagenesis and tectonic deformation.
In the case of carbonate rocks, which are highly soluble in water, additional factors become influential. The flowing groundwater dissolves the flow medium (carbonate rock) Published in the special issue "Five decades of advances in karst hydrogeology" around the existing voids, thus increasing void aperture and the overall hydraulic conductivity of the aquifer. The enlargement of fractures depends on parameters which determine dissolution kinetics (Dreybrodt and Gabrovsek 2003), but mainly the direction and magnitude of groundwater flux (Kiraly 2003).
Considering that groundwater flux depends on hydraulic conductivity and hydraulic gradient, and that hydraulic conductivity depends on the aperture of fractures, the karstification process, which increases void sizes represents an important feed-back loop (autoregulation) in the evolution of karst systems (Kiraly 2003).
As a result, in karst aquifers the hydraulic conductivity field is not determined by the geological history of the rocks alone, but also depends on the dynamic evolution of the groundwater flow system (Kiraly 2003). Another important aspect of the self-regulating process of karstification is that the conductivity field is not constant, but undergoes a dynamic temporal evolution. This evolution is manifested by the hydrodynamic behaviour of karst systems which might indicate various stages of karstification or simply the ability of a system to karstify. Dreybrodt and Gabrovsek (2003) demonstrated, through the use of numerical karst evolution models, that groundwater flux increases slowly at the beginning of karstification, but then is enhanced dramatically when breakthrough of karst conduits occurs. Breakthrough time depends on the initial fracture characteristics (aperture), boundary conditions, and the chemical parameters which determine dissolution kinetics.
The most significant consequence of limestone dissolution associated with karst evolution is increasing hydrogeological heterogeneity. Strong heterogeneity manifests in duality of fundamental hydraulic processes occurring in karst aquifers (Kiraly 1994). Karst duality involves the coexistence of diffuse and concentrated infiltration, flow and discharge.
To develop appropriate investigation techniques and achieve a realistic interpretation of the hydrogeological functioning of a carbonate hydrogeological system, it is necessary to understand the actual degree of karstification and the corresponding hydrodynamic functioning of a system. This requires a process-based quantitative definition and classification of carbonate aquifers. The goal of the present study is to provide a physics-based quantitative classification of carbonate aquifers using hydrograph data, and to demonstrate the applicability of this classification on selected test sites.

Precedents
The definitions and classifications of karst show a wide variety of classification criteria and goals. While most definitions and classifications are based on geomorphological characteristics, there is a very limited number of studies which aim at defining and classifying carbonate systems based on hydrodynamic behaviour.
According to the review of 21 definitions of karst undertaken by Klimchouk (2015), the most important common features included in most definitions are the presence of a specific landscape (karstic landforms), dissolution, and carbonate rocks.
Although these features are typical of karst aquifers, their existence alone is not limited to karstic hydrodynamic behaviour, while on the other hand, it can be present without evidence for some of these criteria. The existing definitions of karst involve several shortcomings: 1. Diagnostic landscape features are not always present-for example, hypokarst can develop without any visible surface manifestations. 2. Dissolution is not limited to karst, because solubility depends on rock type, solute and solvent, pressure, temperature, and pH conditions, and is a widespread phenomenon in nature. 3. Karstic hydrodynamic behaviour is not exclusive to carbonates-for example, evaporite rocks can exhibit karstic groundwater flow (Klimchouk et al. 1996;Gutiérrez et al. 2008;Ford and Williams 2007). 4. Not every carbonate rock produces typical karst landscape features. Dolomites and chalks for example represent ambiguous groups of carbonate aquifers. While some of these aquifers show typical karstic hydrodynamic behaviour (high flow velocities and erratic spring discharges), others completely lack these specific hydraulic characteristics (Sun 1995;Maurice et al. 2006;Mahler et al. 2008;Fournier et al. 2006). 5. Many karst definitions are qualitative and difficult to apply in practice.
Most geomorphic karst classifications have limited hydrologic application because they are mainly based on geological or morphological conditions and most of them have limitations similar to the existing definitions of karst, as they miss the quantitative definition of conditions determining different hydraulic behaviours. Worthington et al. (2017) Worthington et al. (2017) found that none of these definitions are universal, and that definition used may vary with the goal of the study. Some definitions can be useful for characterizing contaminant transport, while others are useful for assessing karst hazards. Besides the fact that the reviewed definitions are qualitative, they are often contradictory or insufficient for a general distinction of karst aquifers.
Classifications based on hydrodynamic behaviour Mangin (1975) introduced a semi-quantitative classification of karst systems based on hydrodynamic behaviour. He applied two parameters for classifying karst systems. The regulation power (K) is defined as a ratio between dynamic volume (the volume of water in storage in the saturated zone above the level of the outflow spring) and spring discharge over a hydrological cycle. Infiltration delay (i) represents the time lag between recharge and spring discharge. Systems exhibiting low i value are characterized by fast infiltration, while systems with i value tending towards one, are recharged by slow or delayed infiltration. The mathematical definition of these parameters is provided in several studies (Mangin 1975;Jeannin and Sauter 1998;El-Hakim and Bakalowicz 2007). Based on these two parameters Mangin (1975) defined five classes as indicated in Fig. 1.
The analysis of spring discharge data from the Milandrine karst system in the Jura, Switzerland, undertaken by Grasso and Jeannin (1998), indicated that although K is stable, parameter i depends on recharge type and catchment size, and thus concluded that the classification of Mangin is uncertain. Later studies by El-Hakim and Bakalowicz (2007) and Liu et al. (2009) also criticized the classification of Mangin (1975) claiming that there are no theoretical reasons for class limits. Hobbs and Smart (1986) presented a qualitative classification and described the main attributes governing the behaviour of carbonate aquifers. Their classification is based on three attributes including recharge, storage and transmission. They claim that these characteristics are independent and they can take different values between end members. In the case of recharge, concentrated and dispersed end members are defined. In the case of flow, conduit flow and diffuse flow are considered as end members. Hobbs and Smart (1986) claim that conduit flow occurs in open channels with high velocities giving rise to turbulent flow. Diffuse flow is confined to tight fractures and pores; velocities are low and flow is laminar obeying Darcy's law. Diffuse flow can occur as intergranular flow in porous limestone or fracture flow in dense limestone. The latter is similar to fracture flow in igneous rocks. The end members describing storage are unsaturated and saturated storage, where the majority of unsaturated storage occurs in the epikarst. The model of Hobbs and Smart (1986) provides a conceptual background for carbonate aquifer classification; however, their model is only qualitative and class limits are not based on actual field measurements. Kiraly (1994) qualitatively defined mature karst as a group of karst systems displaying dual hydrogeological behaviour. According to karstic duality, infiltration, flow and discharge occur both in concentrated and diffuse form. Systems with dual behaviour may be schematized by a high permeability channel network with kilometre-scale meshes, which are immersed in a low permeability fractured limestone matrix, and are discharging through a karst spring. The duality of karst aquifers is a direct consequence of this structure, and thus can be an indicator of karstic hydrodynamic behaviour according to Kiraly (1994).
With respect to the classification of carbonate aquifers and the definition of karst, the most crucial question is how we can quantify the conditions which determine the peculiar hydrodynamic behaviour of karst?
As discussed in the preceding, existing definitions and classifications are mainly qualitative and focus on geological and geomorphological aspects, and most of them lack a quantitative definition for the hydrodynamic classification of carbonate hydrogeological systems.
Although spring hydrographs can provide integral information about the hydraulic behaviour of conduits and matrix blocks, well hydrographs provide information about the dimensions and the hydraulic properties of individual matrix blocks bordered by karst conduits or aquifer boundaries. Although spring hydrograph analysis has an extensive literature, there are only a few studies which discussed well hydrograph analysis for aquifer characterization (Shevenell 1996;Powers and Shevenell 2000;Kovács and Perrochet 2014;Kovács et al. 2015). Jeannin and Sauter (1998) distinguished two groups of spring hydrograph analyses. The same classification of methods can be applied to well hydrographs. Firstly, time series analyses investigate the hydrograph of karst systems during a succession of recharge events. These methods are based on mathematical operations, and their results cannot be directly related to physical characteristics. Although the interpretation of stochastic models can sometimes be related to some physical aspects of the functioning of a hydrogeological system, these methods cannot provide quantitative information about hydrodynamic behaviour.
Single-event methods deal with the hydraulic response of a spring or monitoring well applied to a single recharge event, or with the integrated hydraulic response represented by master recession curves (MRC). An MRC or envelope is a single recession curve derived from a series of recession segments over a wider flow range, and thus providing an average characterization of baseflow response (Kresic 2007).
Most of these methods involve reservoir models or analytical solutions, and are based on physical processes. Analytical methods can be applied for obtaining quantitative information about hydraulic and geometric properties of karst aquifers, and for establishing mathematical relationships between hydrograph features and aquifer properties. Single-event hydrograph analytical methods are applied in this study for establishing a quantitative classification of carbonate hydrogeological systems.

Single-event hydrograph analysis
The plots of spring discharge versus time are referred to as spring hydrographs. The plots of hydraulic head versus time are referred to as well hydrographs. Hydrographs consist of a succession of individual peaks, each of which represents the response of the aquifer given to a recharge event (Fig. 2). Each hydrograph peak consists of a rising and a falling limb. The inflection point on the falling limb indicates the end of direct recharge. The steep segment before the inflection point is called the flood recession, and the flatter segment after the inflection point is referred to as baseflow recession. Quickflow originates from direct recharge through sinkholes and vertical shafts collecting water from sinking streams and from the epikarst. Baseflow consists of water released from the low-permeability matrix blocks, and is much less influenced by the temporal and spatial variations of rainfall than quickflow. The recharge transmitted by matrix blocks takes place through diffuse infiltration from the surface and through gradient inversion between the blocks and neighbouring conduits (bank storage; Kiraly 2003).
Spring hydrographs are traditionally simulated using reservoir models such as those of Mero (1963Mero ( , 1969, Mero and Gilboa (1974) Guilbot (1975), Bezes (1976), Thiéry (2014). The simplest model configuration includes two reservoirs, the first reservoir representing the low-permeability matrix discharging into the conduit system, while the other reservoir representing the high-permeability conduit system. The recent work of Bailly-Comte et al. (2010) pointed out that a parallel configuration of reservoirs is more appropriate than the traditional serial model as it is capable of simulating gradient inversion and the reversal of flow between conduits and the rock matrix (Fig. 3.) The first mathematical models of baseflow recession were provided by Boussinesq (1904) as a hyperbolic form, and by Maillet (1905) as an exponential function. The model of Maillet (1905) is more widely used: where Q t is the discharge [L 3 T −1 ] at time t, and Q 0 is the initial discharge [L 3 T −1 ] at an earlier time, α is the recession coefficient [T −1 ] usually expressed in days. Plotted on a semilogarithmic graph, this function is represented as a straight line with the slope α. This equation is usually adequate for describing karst systems at low water stages. The earlier stages of baseflow recession usually deviate from the simple exponential formulae of Maillet (1905). Forkasiewitz and Paloc (1967) realized that decreasing limbs of hydrograph peaks can usually be decomposed into several (usually three) exponential segments. Forkasiewitz and Paloc During decomposition, the exponential terms are fitted onto the flattest section of a spring hydrograph, and are successively subtracted from the residual hydrograph (Fig. 4). Forkasiewitz and Paloc (1967) assumed that, different segments of a spring hydrograph peak represent different parallel reservoirs, all contributing to the discharge of the spring.
Schoeller (1965) assumed that successive exponentially decreasing limbs on spring hydrographs are due to changes of flow regimes (from turbulent to laminar). Bonacci (1993) discussed that break points on a recession line correspond to changes in a characteristic of the aquifer, such as the decrease of recharge area or the decrease of effective porosity of the aquifer. Kresic (2007) found that recession curve can be approximated by two or three corresponding exponential functions with different coefficients of discharge. He claimed that different exponentials correspond to different microregimes of discharge during recession indicating three types of storage (effective porosity). Kovács (2003) and Kovács et al. (2005) provided an exact analytical solution to baseflow, and Kovács and Perrochet (2008) showed that baseflow recession can be described as a sum of infinite number of exponential components: These authors demonstrated that exponential components do not necessarily represent different permeability classes in an aquifer, but instead they are the consequence of transient phenomena during the emptying of low permeability matrix blocks.

Analytical solutions
The one-dimensional analytical solutions introduced by Boussinesq (1904), Rorabaugh (1964), Berkaloff (1967) and the two-dimensional analytical solutions provided by Kovács (2003) and Kovács et al. (2005) showed that the recession  Forkasiewitz and Paloc (1967) coefficient (α) is proportional to the ratio T /SL 2 where L is block size [L], T is the hydraulic transmissivity [L 2 T −1 ] and S is the storativity [−] of the rock matrix. The work of these authors made the estimation of conduit spacing and rock matrix hydraulic properties possible from hydrograph analysis.
The analytical solutions described in Kovács (2003), Kovács et al. (2005), Kovács and Perrochet (2008), and Kovács and Perrochet (2014) were based on a conceptual model of karst visualized in Fig. 5. This model comprises a rectangular conduit network immersed in the lowpermeability rock matrix forming matrix blocks. The conduit network has one outlet which is the karst spring. This simple but effective model can be described by defining the hydraulic parameters of the low-permeability matrix, the hydraulic parameters of the conduit system, conduit spacing, and the spatial extent of the aquifer (Fig. 5).

Spring hydrographs
An analytical solution for diffusive flux from a twodimensional (2D) asymmetric rectangular block encircled by uniform head boundary conditions can be expressed as follows (Kovács and Perrochet 2008): where initial conditions comprise uniform hydraulic head distribution over the block surface, and where T [L 2 T −1 ] is block transmissivity, S [−] is block storativity, Lx and Ly [L] are block size, β [−] is asymmetry factor. It follows from Eq. (4) that Kovács and Perrochet (2008) demonstrated, that although baseflow can be expressed as the sum of an infinite number of exponential components, higher order components vanish at early times of the recession process, and thus baseflow recession can be sufficiently approximated as the sum of three exponential components (Fig. 6).
The recession coefficients of the first three exponential components of baseflow from an asymmetric block can be expressed as follows: are block dimensions, x and y are distance of observation well from block centre, and α H and α Q are well and spring hydrograph recession coefficients, respectively. After Kovács and Perrochet (2014) For square-shaped blocks, the analytical solution simplifies and the recession coefficients of the first three exponential components of baseflow from a square-shaped block can be expressed as follows: Well hydrographs A one-dimensional (1D) analytical solution for waterlevel changes in a homogeneous block during recession was provided by Rorabough (1960), who aimed at using it for the estimation of aquifer diffusivity from well hydrograph data. Besides the work of Rorabough (1960), very few studies discuss the application of well hydrographs for parameter estimation in karst. These include Atkinson (1977), Shevenell (1996) and Powers and Shevenell (2000). These studies, based on 1D solutions, were mainly focussed on the estimation of hydraulic parameters, and did not attempt to estimate block geometry of karst aquifers from well hydrograph data. Kovács and Perrochet (2014) provided a 2D analytical solution for hydraulic head evolution in a rectangular matrix block during baseflow, and tested and demonstrated the solution on synthetic numerical models. A field application of the 2D analytical solution of Kovács and Perrochet (2014) for the estimation of hydraulic properties and conduit spacing was undertaken by Kovács et al. (2015).
The 2D analytical solution of Kovács and Perrochet (2014) was based on the conceptual model indicated in Fig. 5. The process of baseflow recession in matrix blocks surrounded by karst conduits is indicated in Fig. 7.
The 2D analytical solution of Kovács and Perrochet (2014) can be expressed as: It can be seen from the analytical solution, that the falling limb of the hydrograph of a well installed in a lowpermeability matrix block can be expressed as the sum of an infinite number of exponential terms similarly to spring hydrograph baseflow. The recession coefficients of the first three exponential components of well hydrograph baseflow can be expressed similarly to spring hydrographs making use of Eqs. (7), (8) and (9). As a consequence, the baseflow decomposition concept introduced by Forkasiewitz and Paloc (1967) can be equally applied to well hydrographs. The details of well hydrograph decomposition are discussed in Kovács and Perrochet (2014).

Quantitative classification
A sensitivity analysis using numerical combined discrete-continuum models (Király 1979(Király , 1985Kiraly and Morel 1976a) was undertaken by Kovács (2003) and Kovács et al. (2005) in order to investigate the influence of aquifer parameters on the hydraulic behaviour of karst systems. This study also established a quantitative relationship between hydraulic parameters and recession coefficient, thus making the estimation of aquifer properties from hydrograph data possible. The analytical solutions and sensitivity analyses introduced by Kovács (2003) and Kovács et al. (2005) Fig. 6 Exponential components of the analytical solution of Kovács and Perrochet (2008) established the basis for the development of a quantitat i v e c l a s s i f i c a t i o n o f s t r o n g l y h e t e r o g e n e o u s hydrogeological systems including carbonate aquifers.
The sensitivity of the recession coefficient of the flattest section of baseflow (α 1 in the analytical models) to aquifer hydraulic properties including conduit conductivity and conduit spacing were investigated by Kovács (2003) using the conceptual model indicated in Fig. 5. During the sensitivity analysis square-shaped matrix blocks were applied in order to reduce the number of parameters. Conduit storage coefficients S c [L] were calculated from relevant conduit apertures assuming water compression. Figure 8 indicates that increasing the value of conduit conductance generates increasing values of the recession coefficient (α 1 ) until it reaches the value expressed by the analytical solution for block recession (Eq. 10). This is the highest value of baseflow recession coefficient α 1 achievable by the system during karstification assuming constant hydrogeological conditions. Conduit diameter can be expressed from the applied values of conduit conductance by applying Poiseuille, Strickler or Louis formulas. Conduit diameters smaller than the threshold value throttle the flow and hinder the free emptying of matrix blocks, thus resulting in lower recession coefficients. Among these conditions, the entire system (catchment) exhibits diffuse or in other words divergent flow. This situation is valid during the early stages of karstification, but it is also an intrinsic characteristic of systems which are not sufficiently karstifiable. Because of the influence of conduit/fracture aperture on the recession process, this flow condition was defined as conduit-influenced flow regime (CIFR) by Kovács (2003). In this study, the CIFR condition is simply called divergent flow, following a standard terminology. A mathematical expression of recession coefficient under CIFR conditions was provided by Kovács (2003) and Kovács et al. (2005) as: , T m is matrix transmissivity, S m is matrix storativity, f is conduit frequency (f = 1/L where L is block size of a square-shaped block). After exceeding the threshold value, the increase of conduit conductance has no further influence on the recession coefficient; the recession process is restricted by the hydraulic and geometric parameters of the low-permeability matrix blocks alone, and Eq. (4) provides an adequate characterization of the systems global response. At this stage the hydraulic gradient in the conduits becomes negligible, and conduit flow has no further influence on the drainage of the low-permeability blocks. The conduits act as constant-head boundary conditions as assumed by the analytical models (Eqs. 4 and 13). This flow condition was defined as matrix restricted flow regime (MRFR) by Kovács (2003). In this study, the MRFR condition is simply referred to as convergent flow.
The convergent flow condition is an evident consequence of the positive feedback loop of the karstification process. Conduit diameters increase until they reach sufficient conductance for draining the water discharged from the low-permeability matrix blocks. The conduit system of a mature karst which has reached hydrodynamic equilibrium through conduit enlargement is capable of draining the entire baseflow discharge received from the rock matrix. This situation usually entails open-channel flow in the conduits. On the other hand, the same system may become flooded and the drainage capacity of conduits exceeded during flood events (Bailly-Comte et al. 2010). This situation entails fully constricted channel flow during flood conditions. The condition of gradient inversion is discussed in (Kiraly and Morel 1976a, b;Kiraly 2003). This also means that turbulent vs. laminar flow has influence on the temporal variations of spring discharge only during floods. Conduit flow type has no influence on the recession coefficient after the end of flood discharge in the case of convergent flow; it is controlled by the drainage of matrix blocks alone.
The conceptual characterization of the connection between karstification and recession characteristics can be summarized as follows. The gradual increase of conduit diameter is supposed to happen along the curve depicted in Fig. 8a and outlines the evolution of the recession coefficient during the karstification process. As discussed by Dreybrodt and Gabrovsek (2003), karstification evolves slowly during the dissolution process. During this stage the diameter of conduits slowly increases, and there are no karstified conduits connecting the recharge area to the karst spring. This situation corresponds to the divergent or CIFR flow condition defined by Kovács (2003). At later times, an abrupt increase in discharge rate happens, when breakthrough of the first karst conduit occurs (Worthington and Ford 2009). After breakthrough happens, the karstification process speeds up, and corresponding spring discharge increases. Besides the master conduit, transverse conduits develop, which decrease the representative matrix block size, and further increase the recession coefficient. The process continues until the system reaches equilibrium, and the final form of the conduit network is established. At this stage the representative block size is reached, and the entire system manifests convergent or MRFR type hydraulic behaviour during baseflow.
It is demonstrated in the preceding how the karstification process influences the values of the recession coefficient. It is also shown that two significantly different flow conditions may exist in carbonate rocks (and in strongly heterogeneous systems in general).
In systems with well-developed conduit networks, the baseflow is controlled by the drainage of matrix blocks. This flow condition is referred to as convergent flow. In these systems the baseflow recession is influenced by the size and the hydraulic properties of the low-permeability matrix blocks alone. The dominant flow direction in the rock matrix during baseflow is towards karst conduits, which collect diffuse discharge from matrix blocks and transfer it towards the karst spring. This is considered to be a crucial criterion of karstic hydrodynamic behaviour and an ultimate consequence of the karstification process. These systems are also considered to satisfy the qualitative definition of duality introduced by Kiraly (2003), as only well-developed conduit systems can receive, transfer and release concentrated flow.
In systems with well-developed karst conduits, baseflow is only influenced by the hydraulic properties and the geometry of matrix blocks. This convergent flow condition is defined here as a typical characteristic of well-developed karst systems manifesting karstic hydrodynamic behaviour.
In systems with poorly developed karst conduits or unkarstified fractures, the baseflow recession is influenced by the aperture size and hydraulic parameters of fractures/ conduits. Also, the hydraulic properties of low-permeability blocks, fracture spacing, and aquifer extent have influence on baseflow recession. This divergent flow condition is considered to be a characteristic behaviour of poorly karstified systems and fissured systems, manifesting nonkarstic hydrodynamic behaviour.
The fundamental difference between these two types of carbonate hydrogeological systems is indicated in Fig. 9. Whereas both systems in the figure have identical matrix hydraulic properties and identical conduit spacing, the first system has higher conduit conductance than the second system.
From the hydrodynamic point of view, hydrogeological systems manifesting convergent baseflow can be defined as karst systems, whereas systems manifesting divergent baseflow can be defined as fissured systems. This classification is irrespective of geomorphological or geological conditions, and is purely based on hydrodynamic behaviour.

Classification methodology
The classification of a carbonate system which is based on hydraulic behaviour requires a systematic analysis of hydrograph data and careful comparison with field observations. In order to   Kovács and Perrochet (2014) classify carbonate systems, all of the available data need to be analysed, including information about conduit spacing based on surface manifestations, cave maps, geophysical data, and hydraulic properties of the rock matrix obtained through pump tests. Hydrograph analyses of spring and well hydrographs, whichever are available need to be undertaken. The simultaneous analysis of spring and well hydrographs provides the most valuable information about a carbonate aquifer. The classification may be made on the basis of estimation of geometric or hydraulic properties from hydrograph data, by making an assumption about the dominant flow regime (convergent or divergent) and comparison of the results with field observations. A satisfactory accordance between calculated properties and field observations may confirm the initial assumption about flow regime, while a significant mismatch will disprove it. This method supports a reliable hydrodynamic classification of a carbonate system and the selection of proper investigation methods and interpretation of flow conditions.
Parameter estimation can be achieved through the following steps: 1. Preparation of a master recession curve or selection of hydrograph sections with long undisturbed recession periods. 2. Hydrograph analysis and identification of α 1 , α 2 , and if possible α 3 baseflow recession coefficients by curve fitting. 3. Calculation of the asymmetry factor β from α 1 and α 2 : Calculation of block size L x and L y from α 1 and β: Or, calculation of diffusivity D = T/S from block size and recession parameters:

Field examples
In order to demonstrate the applicability of the quantitative classification introduced in the previous section, the proposed methodology was applied to two different carbonate systems. These systems represent typical examples for carbonate aquifers which are utilized as drinking water resources, and thus require a detailed understanding of flow dynamics.

Szinva catchment, Bükk Mts., Hungary
The Szinva catchment extends over the Fennsík (uplands) zone of the Bükk Mountains and covers an elongated area of 12 × 3 km (Fig. 10). Catchment delineation was undertaken by Hernádi et al. (2013) based on tracing tests, surface topography, and water budgets. The upland area geology comprised middleand upper-Triassic limestones. The main area of the Szinva catchment is mainly hosted by the anchimetamorphosed platform facies limestones of the Bükkfennsík Formation. The south-eastern part of the catchment is composed of the Felsőtárkányi Formation, which is thick bedded, intrashelf basin facies limestone including marl and chert intercalations (Miklós et al. 2020; Fig. 10).
The uplands area has a large number of sinkholes (Hernádi et al. 2013) which follow linear features that indicate the location of karst conduits at depth, as discussed by Kovács et al. (2015). The horizontal distance between suspected conduits outlined by sinkholes is between 200 and 600 m (Fig. 11).
The outlet of the Szinva catchment is the Szinva spring, which is located at an elevation of 340 masl. Karst water levels are observed within the catchment at two observation wells. Observation borehole NV-17 is located 8,500 m west of the Szinva main spring and has an average water level of 525 masl. This borehole is screened in Triassic age Bükkfennsík Limestone. Observation borehole M-6 is located 4,900 m west of the Szinva spring and has an average water level at 450 masl. This bore is also screened in Triassic age Bükkfennsík Limestone.

Hydrograph analysis
The investigation of the Szinva catchment included the analysis of well hydrographs of monitoring boreholes Nv-17 and M-6 and also the analysis of the Szinva spring hydrograph. The aim of the analysis was to derive information about the representative block size on the basis of hydraulic behaviour and its comparison with field observations about suspected conduit network geometry to confirm or disprove assumptions about aquifer classification.
Individual hydrograph peaks followed by long recession periods were selected for analysis. Master recession curves generated from entire time series were also analysed. Monitoring bore NV-17 has an hourly dataset of 22 years ranging between 10/10/1992 to date, while monitoring bore M-6 has a short dataset comprising 21 months of daily data. Some sporadic water-level data were also measured between 1995 and 2001. The excessive sharp peaks prevailing on the M-6 hydrograph are believed to indicate the influence of a karst conduit intersecting the bore profile. The Szinva spring has a 21-year-long water-level time series. Water levels were measured in an artificial spring basin with a prismatic geometry. Hydrograph recession in such settings has similar characteristics and can be described with the same recession coefficients as spring discharge. For this reason, in the case of the Szinva spring, water level time series were analysed. The hydrographs of monitoring bores NV-17 and M-6 and of the Szinva Spring are shown in Fig. 12.
In this study a manual approach to hydrograph analysis was applied, which included the manual fitting of the flat exponential baseflow segment on semi-log head vs. time plots, and the subtraction of the fitted curve from the hydrograph as suggested by Forkasiewitz and Paloc (1967). The process was repeated until a reasonable fitting was possible. In most cases two or three exponentials could be fitted on the hydrographs.
The recession coefficients obtained from the decomposition of selected peaks and master recession curves of these hydrographs are summarized in Table 1. The analysis was undertaken according to the formulae introduced in the previous sections.

Discussion
Hydraulic conductivity of the rock matrix was determined from well tests. An average value of K = 1 · 10 −5 m/s (Székely et al. 2019) was applied in the analysis. The specific yield of the low-permeability rock matrix was determined from well tests and literature data, and an average value of S y = 0.005 was applied.
On the bases of hydrograph decomposition, three ranges of recession coefficient were identified:4 · 10 −11 < α 1 < 7 · 10 −10 L/s, 1 · 10 −7 < α 2 < 1 · 10 −6 L/s, and 1 · 10 −6 < α 3 < 2 · 10 −5 L/s. It can be seen in the preceding, that α 1 is two or three orders of magnitude lower than α 2 . According to the analytical solutions assuming realistic block shapes, the ratio between α 1 and α 2 should remain within one order of magnitude. As a conclusion, the significant difference between α 1 and α 2 in the Szinva catchment cannot be explained by the recession of the same aquifer volume.
The calculation of block size from the late recession coefficient α 1 indicates the drainage of an extensive aquifer zone (L = 22-31 km). NV-17 is located on the assumed catchment boundary, and on the basis of this analysis it indicates the recession of a smaller, hydraulically isolated catchment of L = 7-8-km size. The calculation of geometry from intermediate (α 2 ) and early (α 3 ) recession coefficients indicate block sizes between L = 200-600 m.
Hydrograph analysis reveals the temporal and spatial changes in the hydraulic functioning of the investigated karstic catchment, and contributes to the conceptual understanding of its hydrogeology. Results of hydrograph analysis indicates, that after long periods of drought, water table drops below the level of active karst conduits, and block-scale drainage of the rock matrix is replaced by regional nonkarstic flow (Fig. 13). These dynamic changes in flow scale over the recession period are expressed in the variations of recession coefficient of spring and well hydrographs Kovács et al. (2015).
The domain size of 24-31 km calculated from the late recession coefficient α 1 corresponds well with the actual size of the Bükk carbonate plateau (~20 km). Calculated block size (L~200-600 m) is in accordance with expectations based on field observations, which confirms the assumption about the dominant convergent flow domain. It is concluded from this study that the Szinva catchment behaves hydraulically as a mature karst system.

Kádárta catchment
The village of Kádárta belongs to the city of Veszprém, in western Hungary (Fig. 14). The waterworks of Kádárta supplied water to Veszprém from two spring galleries located at Fig. 9 Classification of carbonate systems based on hydrodynamic behaviour. Blue lines indicate equipotentials, while small red arrows indicate flux vectors. Classification is determined on the basis of hydrodynamic behaviour. Premature karsts are classified as fissured systems. T m [L 2 T −1 ] is matrix transmissivity, S m [−] is matrix storativity, α 1 is baseflow recession coefficient [T −1 ], K c is conduit conductance [L 3 T −1 ], f is conduit frequency (f = 1/L where L is block size), and A [L 2 ] is catchment area the northern edge of a carbonate plateau made up of Triassic platform dolomite. Hydrogeological investigations were performed by Kovács (1998), Csepregi et al. (1998 and Kovács et al. (2001) following nitrate contamination identified in groundwater, which was concluded to originate from agricultural fertilisers applied over the recharge area.
The Veszprem plateau consists of NE-SW tending strips of upper-Triassic calcareous sediments, dipping 15-30°to the NW (Fig. 14). Thrust faulting along two major reverse faults resulted in the recurrence of these strips in the area. The northern fault line, named the Veszprém thrust zone, represents a hydraulic barrier and the north-western boundary of the Kádráta aquifer.
The Kádárta aquifer consists of white, well stratified, middle-Triassic Budaörs dolomite (Budai 2006). The aquifer is hydraulically bordered by the outcrop of underlying mid-Triassic low-permeability calcareous sediments on the SE, and the outcrop of overlying upper-Triassic marls from the NW. The dolomite aquifer pinches out on the SW in the area of Veszprémfajsz.
The hydraulic gradient is about 10-12 m/km, with a flow direction to the N-NE. The Kádárta aquifer discharges at several springs, of which the Kádárta springs have the highest discharge. These springs are located at the intersection of a N-S trending traverse fault and the Veszprém thrust zone. The surface manifestation of this transverse fault is a dry valley called Kasza Valley.
The Kadarta-springs consist of two springs (western and eastern spring), located on both sides of the bottom of the Kasza Valley. The total discharge of the Kádárta springs is about 8-11,000 m 3 /day. A numerical groundwater flow model by Kovács (1998) outlined the recharge area of the springs and estimated an average groundwater travel time of 20-30 years in the Veszprém plateau area.   Table 1 Results of hydrograph decomposition, Szinva catchment, Bükk Mountains, Hungary. α 1 and α 2 and α 3 refer to the recession coefficients in Eq. (3). L block and L catchment values were estimated using the symmetrical analytical solution (Eqs. 10-12). L x and L y were estimated using the asymmetric analytical solution (Eqs. 7-9). Average values in italic The Kádárta dolomite aquifer is approximately 10 km long, 2 km wide, and its estimated thickness is up to 1,000 m. The dolomite is covered by Quaternary loess in the SW area of the catchment. The maximum thickness of loess is 5-7 m, and it is absent in the northern part of the catchment. The soil is thin (10-20 cm) where it overlies the dolomite surface.
A surface geophysical study to investigate fracture density was undertaken by Kuslits (2014). This study used various geophysical methods including radio-magnetotellurics (RMT), ground penetrating radar (GPR) and electrical resistivity tomography (ERT). The study concluded that the fracture spacing of NW-SE trending faults on the plateau is between 150 and 200 m in the SW-SE direction. No significant SW-NE trending fractures were identified besides the major tectonic elements of the area such as the Veszprém Thrust.
Additional geophysical surveys were undertaken by Szalai et al. (2018) to compare various geophysical methods such as ERT and pricking-probe (PriP). A detailed survey was undertaken in the vicinity of the Kádárta springs. A number of individual fractures were localized with a spacing between 1.5 and 15 m.

Hydrograph analysis
The monitoring network comprises six monitoring bores. Within the framework of this study, three of monitoring bores were equipped with automatic pressure transducers. Five-yearlong hourly well hydrograph data were collected. The boreholes are 25-35 m deep with screen intervals between 20 and 30 m. The hydrographs of monitoring boreholes F-1, F-2 and F-3 were analysed with the aim of deriving characteristic fracture/ conduit spacing. Pumping tests were undertaken in each monitoring borehole to determine hydraulic properties. The average hydraulic conductivity determined from pump tests is K = 1.7 · 10 −5 m/s. A specific yield of 0.08 was used in this study.
To measure discharge of the spring, a weir was installed in the spring gallery. Spring discharge was measured from May 2015 to October 2018. There is a data gap between June 2016 to May 2017 because of instrument failure. Discharge was calculated from the appropriate rating curve.
The hydrographs of the monitoring boreholes show a synchronized fluctuation of water levels. Borehole F-2 in the Kasza Valley shows a lower groundwater level. The hydrograph of the Kádárta spring shows significant tidal fluctuations. This is believed to be a consequence of the fissured nature of the aquifer and the tidal variations of fissure aperture and thus of the effective porosity of the aquifer. The long-term variations of the spring hydrograph show a behaviour similar to that of the wells, but the falling limbs of the hydrograph peaks are steeper than those of the well hydrographs. The hydrographs of monitoring boreholes F-1, F-2 and F-3 and of the Kádárta spring are indicated in Fig. 15.
While the mathematical expressions of recession coefficients are identical for spring and well hydrographs (Eqs. 7-9), the decomposition technique applied for well hydrographs may differ slightly from those of the spring hydrographs. For spring hydrographs, the exponential components are always positive, but for well hydrographs, the second component can be either negative or positive depending on the position of the well within a matrix block.
If the second exponential component (H 2 ) is negative (convex hydrograph peak), a complementary decomposition approach can be applied. This includes the fitting of the first exponential on the flattest section of a well hydrograph, and the subtraction of the time series from Fig. 13 a Conceptual hydrograph, b Conceptual flow model of the Szinva catchment, Bükk Mts., Hungary. 1 = during flood recession, when conduits feed the rock matrix, 2 = during early baseflow recession, when the conduits drain the rock matrix, and 3 = during late baseflow recession, when the water level drops below the karstification horizon, and flow-scale increases dramatically this fitted exponential function. A second exponential can then be fitted on the residual hydrograph. The decomposition of selected sections of the well hydrographs of wells F-1 to F-3 and the Kádárta spring hydrograph are demonstrated in Fig. 16.
The well hydrograph peaks obtained from boreholes F-1, F2 and F-3 could be decomposed into two exponential segments. The hydrograph peaks of the spring gallery could not be decomposed and thus only one fitted exponential curve is indicated in Fig. 16d. The results of hydrograph decomposition are summarized in Table 2.

Discussion
It can be seen in Table 2 that the recession coefficient α 1 of the analysed spring hydrograph is two orders of magnitude higher than those of the well hydrographs. As only one exponential component could be identified on the spring hydrograph, the expression of recession coefficients based on the symmetric block assumption was applied for comparison between the results.
The representative block size estimated from the spring hydrograph peaks is approximately 500-700 m. The block size estimated from the well hydrographs is approximately 6-9 km. This indicates an order of magnitude difference in estimated block size when using different data sources (well hydrographs vs. spring hydrographs).
The significant difference in block size estimation might have the following reasons: 1. The spring is draining an aquifer volume with significantly lower block size than that indicated by the monitoring wells. This is a realistic assumption considering that the aquifer is subdivided into two blocks by the transverse fault manifested by the Kasza Valley dissecting the dolomite plateau. 2. The spring is draining an aquifer volume with significantly higher hydraulic diffusivity than that indicated by the monitoring wells. This is possible if we assume a higher fracture density around the spring. The results of surface geophysical survey support this assumption. This can be the result of the tectonic setting: The spring is located at the intersection of the Veszprém thrust and the Kasza Valley transverse fault. A subzone within the aquifer with higher diffusivity also means that the size of this zone is smaller than that of the entire aquifer. This assumption implies two blocks with different sizes and hydraulic diffusivity values.
In order to test the aforementioned scenarios, an analytical modelling experiment was undertaken. The setup of the model is indicated in Fig. 17. The symmetrical block formulation was applied for the calculation of disc h a r g e b y s u m m i n g s i x e x p o n e n t i a l b a s e f l o w components. The assumptions about block size and hydraulic properties together with corresponding recession coefficients are listed in Table 3.
Model results suggest that the aquifer is in fact subdivided into two blocks by the transverse fault manifested by the Kasza Valley dissecting the dolomite plateau (Fig. 14). While the larger western block includes the monitoring boreholes, the small eastern block has no monitoring boreholes installed in it. As scenario 2 (in Table 3) provides a better fit to observations, it can be assumed that the small block is also more intensely fractured because of the vicinity of large structural elements, and thus has higher hydraulic diffusivity.
The Kasza Valley behaves as a boundary condition for flow discharging from the two dolomite blocks located east and west of Kasza Valley. This is evidenced by the water-table topography and the presence of temporary springs along the Kasza Valley bottom. The Kádárta springs discharge water originating from both blocks.  Table 2 Decomposition results from the Kádárta aquifer. α 1 and α 2 refer to first and second recession coefficients in Eqs. (3). L block values were estimated using the symmetrical analytical solution (Eqs. 10-12). L x and L y were estimated using the asymmetric analytical solution (Eqs. 7-9) The conceptual model of the Kádárta dolomite aquifer is indicated in Fig. 17. According to this model, on the bases of hydrograph analysis, the aquifer contains two matrix blocks of different size and different hydraulic properties. The aquifer is drained by the Kádárta springs, and the larger block contains monitoring wells which provide information about its hydrodynamic behaviour. Spring discharge can be represented by combining the base-flow discharges from the two blocks. Inverse analytical modelling suggests that the measured hydrograph parameters can be best explained by the assumption of higher hydraulic diffusivity in the small block, an assumption supported by field observations.
With regard to the classification of the Kádárta aquifer, on the basis of hydrodynamic behaviour, it is evident that the fractures identified through geophysical surveys do not behave as karst conduits. Although the fracture spacing identified in the small block is between 1.5-15 m, the representative block size derived from hydrograph decomposition is between 500 and 1,000 m. The same applies to the larger block-while fracture spacing determined over the plateau from Radiomagneto-Tellurics is between 150 and 200 m, hydrograph analysis yielded a representative block size between 6 and 9 km.
The only feature behaving as a karst conduit in the Kádárta catchment is the Kasza Valley transverse fault. This feature dissects the aquifer into two blocks and behaves as a boundary condition for baseflow. This feature is assumed to be a tectonic feature rather than the result of karstification.

Summary and conclusions
The existing classifications of carbonate hydrogeological systems are mainly based on geological or geomorphological conditions. Karstic groundwater flow is not always adjusted to surface geomorphology, and geology itself does not determine karstic hydrodynamic behaviour. This study introduces a quantitative classification based on hydrodynamic behaviour, which is necessary for understanding the physical functioning of a carbonate hydrogeological system, to predict future conditions, and to select appropriate investigation techniques.
Analytical methods based on single hydrologic events are illustrated in this study for establishing a quantitative classification of carbonate hydrogeological systems. The methodology was developed for both spring and well hydrographs. Baseflow recession coefficient is used as an indicator of hydrodynamic behaviour. Analytical solutions were developed in order to establish a link between recession coefficient and aquifer characteristics which serve as a basis for the classification.
Sensitivity analyses performed on synthetic numerical models identified two different flow regimes, and revealed the transition between fissured and karstic flow dynamics during the process of karstification. Analytical solutions were applied for characterizing these flow conditions and for defining the quantitative criteria for separating the two representative flow regimes. The threshold value of conduit conductance for quantitative classification is: This formula expresses the quantitative classification of carbonate systems.
Conduit conductance smaller than K c * throttle the flow and hinder the free emptying of matrix blocks during baseflow. This situation exists during the early stages of karstification, but it is also an intrinsic characteristic of systems which are not karstifiable.
Conduit systems with a conductance higher than K c * have no limiting influence on baseflow recession coefficient. The dominant flow direction in the rock matrix during baseflow is towards karst conduits. This is considered to be a crucial characteristic and the criterion for a qualitative classification of karstic hydrodynamic behaviour and an ultimate consequence of the karstification process. In these systems the recession process is restricted by the hydraulic and geometric parameters of the lowpermeability matrix blocks alone. This flow condition is a typical indicator of mature karst systems which have reached hydrodynamic equilibrium. The hydrodynamic classification of a carbonate system can be undertaken through the estimation of its hydraulic and geometric parameters. As all of these parameters are rarely available, an alternative way of classifying a carbonate system based on its hydrodynamic behaviour is by making an assumption about the dominant flow regime, estimating the geometric properties of the conduit system from hydrograph data, and comparing the results with field observations. In order to demonstrate the applicability of the quantitative classification, the proposed methodology was applied to two different carbonate systems.
At the Szinva spring catchment located in the Bükk Mountains, NW Hungary, hydrograph analysis identified three ranges of baseflow recession coefficients both on spring and well hydrographs. While α 2 and α 3 indicated block-scale baseflow recession, α 1 indicated the recession of a regional unkarstified aquifer volume. Hydrograph analysis revealed a change of flow scale triggered by the drop of the groundwater table below the karstified zone. The block sizes calculated from spring and well hydrographs agree with suspected conduit spacing derived from sinkhole locations. It is concluded from this study that the Szinva catchment behaves hydraulically as a mature karst system.
At the Kádárta dolomite aquifer located in West Hungary, the analysis of well hydrographs yielded block sizes significantly larger than those calculated from spring hydrographs. Results of hydrograph analysis were compared with the outcomes of surface geophysical surveys. A systematic analysis of geological and hydrogeological field data suggest that the aquifer contains two matrix blocks of different size and different hydraulic properties. These blocks are separated by a tectonic feature dissecting the dolomite aquifer. Inverse analytical modelling confirmed this conceptual model. The fractures identified through geophysical surveys do not behave as karst conduits. The only feature behaving hydraulically as a karst conduit is the Kasza Valley transverse fault.
The appropriate classification of dolomite aquifers is particularly important, and requires caution in the selection of investigation tools and interpretation of hydrogeological data. These systems might exhibit either karstic or fissured type hydrodynamic behaviour which cannot be determined purely based on geology and/or geomorphological analysis. Hydrograph analysis provides a reliable means of determining the hydrodynamic functioning of an aquifer and supports the quantitative classification of carbonate hydrogeological systems.