Advances, challenges and perspective in modelling the functioning of karst systems: a review

We present a discussion of the state-of-the art on modelling geometrical characteristics, hydrogeological behavior and geochemical evolution of karst aquifers of meteoric origin. The considered key karst processes include: (1) the recharge processes, (2) the various hydrodynamic processes in the vadose and phreatic zones, (3) the related transport processes and (4) the speleogenesis processes. Different types of approaches for modelling geometrical characteristics of karst conduit networks are summarized. Integrated numerical studies on hydrogeological behavior of karst aquifers with functional and physically based models are then reviewed. Challenging issues in karst modelling are further discussed based on previous modelling progresses. The paper finally offers recommendations for advancing the modelling of hydrogeological behavior of karst systems and lists several open questions for future research.


Introduction
Karst aquifers are highly heterogeneous media that are composed of numerous flow paths of high permeability formed by the dissolution of carbonate rocks. The primary porosity of the carbonate rocks can play a major role in the karst genesis processes, as well as the fractures and bedding planes within the rock mass, which are often enlarged by dissolution, forming hierarchically organized flow networks.
Modelling key karst processes is thus a quite challenging issue. This difficulty ultimately hinders our ability in water management of karst groundwater resources and prediction, prevention and mitigation of aquifer contamination.
The first challenge in modelling karst is to solve the complex hydrodynamics problem, i.e., the physical processes of fluid flow and transport in karstic rocks. Karst aquifers consist of multiple media that conduct flow (e.g., conduit, fracture and matrix) with distinct hydraulic properties. Accordingly, the hydraulic conductivity may vary by several orders of magnitude over a short distance where it exhibits a sharp contrast and strong anisotropy (Kuniansky et al. 2012). In addition, the various flow features are hierarchically organized such that multi-scale flow behaviors are very common. As a result, the behavior of karst aquifers is characterized by the striking characteristic of duality (Goldscheider and Drew 2014). Although numerical techniques exist for numerical modelling of flow processes in individual compartments, the coupling of different flow processes remains challenging. The second fundamental challenge in modelling karst is the geometrical representation of complex three-dimensional karstic conduit systems. In most cases, karst forms by dissolving limestone fractures, which occur at all scales. Though some large cavities may be detected by geophysical tools, direct observation of the detailed 3D structure/geometry of karst conduit networks in the subsurface is generally not possible. The traditional approach of direct data sampling relies on cave exploration, outcrop cross-section mapping and borehole logging. Data collections by caving and outcrop mapping are time consuming and remain limited. Considering the technological advances in well-logging and in-situ testing, a detailed, continuous, 1 3 396 Page 2 of 26 characterization of the distribution of discontinuous features (e.g., fractures, bedding planes, vugs, karstic conduits, etc.) is possible. However, the resulting data are largely limited to one-dimensional and only within the regions adjacent to the observation boreholes. Seismic surveys may be able to provide spatial dimensions about the 3D large-scale features such as faults and cavities, but their application in karst media is often limited due to a resolution threshold. Consequently, the 3D description of karst geometrics requires extrapolating 1D/2D observations in small samples of the karst reservoir to the whole study domain in 3D. Therefore, the question of how to create realistic karst networks remains an unresolved issue.
Another critical issue of karst modelling concerns the hydrological and hydrochemical functioning of karst media under complex artificial or natural hydraulic conditions. Groundwater movement and contaminant migration through karst systems is predominately governed by conduits and caves. Recharge is often highly focused at sinkholes or stream networks and therefore is generally less dependent on evapotranspiration and soil cover than in other hydrogeologic settings. Discharge from karst aquifers primarily occurs at conduit-flow karst springs, and can be very large and highly variable both spatially and temporarily. The focused discharge from springs, however, represents an integrated response of water movement through different flow paths and mixing processes that occurred in various locations within the karst system. The understanding of hydrochemical changes in karst aquifers is also essential for groundwater protection and aquatic ecosystem management. Given these complexities and challenges, many numerical approaches have been proposed in the last decades to simulate the main karst processes. All methods have both advantages and disadvantages. There is no a systematic and complete numerical model that can accurately describe the flow and transport processes in karst aquifers. Relevant reviews of numerical approaches for modelling karst flow processes exist (Ghasemizadeh et al. 2012;Kovács and Sauter 2007). Their main focus was to classify karst flow models and the various numerical approaches that, in most reviews, were divided into two general classes based on the type of parameterization: the lumped and distributed parameter models. Even though distributed lumped parameters approaches were proposed to account for the spatial heterogeneity of either epikarst (Hartmann et al. , 2015 or land cover land uses (Bittner et al. 2018) at the scale of the catchment's scale, the lumped parameter models generally do not consider spatial variability of the karst system hydrodynamic properties and rather focus on the temporal variation of global precipitation-recharge-discharge behavior of the system. In contrast, the distributed parameter models aim at mimicking the distributed physical properties that controls both the spatial and temporal variations of hydraulic response of the system. This classification was useful for a general category but may become irrelevant for the new generation of models. Indeed, some models now couple a functional approach describing the flow process in one system component and a physically based approach that captures the spatial variability of another component. In addition, novel modelling methods that introduce simple spatial discretization (often one-dimensional variation) into the functional approach also emerges in recent years. For these reasons, in this review, we do not refer to "lumped" and "distributed parameter models" but rather adopt the terminology of 'functional approach' and 'physically based approach' when classifying models.
In previous reviews, the various models were discussed separately, which did not allow for a thorough presentation of the links between them, let alone the possibility for model and data integration. Here, we attempt to provide an integrated point of view on karst modelling and to better describe the links between the various models. Accordingly, we put emphasis on the connections and combinations of existing models for an improved characterization of karst systems, as well as on the use of data from various sources through the state-of-the-art inversion techniques. The karst modelling approaches that will be discussed in this review paper are summarized in Fig. 1. Two major group of models are considered: (1) models for karst hydrology and hydraulics simulation and (2) models for karst generation and geomorphology simulation. Both forward and inverse modelling techniques are reviewed. We highlight the models and techniques that emerge in recent years (those indicated by red font color in Fig. 1). We also discuss the specific issues than each model can address, and their potential impact on the management decision making. This paper focuses on discussing numerical modelling approaches for karst processes simulation. It is structured in the following way. "Key karst processes and related conceptual models" details the key processes occurring in epigene karst systems and highlights the importance of accurately describing them in conceptual models before capturing them in numerical models. In "Modelling of key karst processes" and "Modelling of karst conduits network", we present recent advances in numerical modelling of key karst processes and karst network geometry, respectively. In "Modelling of karst hydrogeological functioning", we show what can be achieved in terms of simulation of karst hydrogeological functioning with the current modelling tools and discuss the limitations of these tools. We finish this review by providing several perspectives for future modeling approaches, with the aim to better understand, characterize and reproduce key flow and transport processes that occur in natural epigenic karst systems.

Key karst processes and related conceptual models
The first step in karst modelling is the building of an appropriate conceptual model (Fig. 2). The conceptual model provides a collective visualization of the hydrodynamic cycle: recharge, flow/storage and discharge. A well-structured conceptual model should be readily translated into a mathematical model formulated with a set of coupled differential equations. There exist several conceptual models for karst aquifers in the literature. A classical conceptual model originates from the dual-medium model developed for fractured petroleum reservoirs (Drogue and Bidaux 1992). In this model, the main flow paths correspond to fracture and karst conduits while the fissured matrix contributes to the main storage. The model was quite successful for representing flow processes in the phreatic zones of karst aquifers but it is hardly applicable for karst systems where the unsaturated zone processes are significant. By examining the spatial and temporal variation of water table level and tracer concentration, it was found that the shallow portion of the karst aquifer, i.e., the soil zone and epikarst, may play an important role in regulating flow and transport processes (Perrin et al. 2003). The concept of epikarst as a shallow, highly karstified zone below the soil zone was introduced by Mangin (1975) and Williams (1983). The epikarst was also described as playing a major role in the storage of karst aquifer (Klimchouk 2000;Sauter 1992;Williams 1983) and in the sequence of vertical infiltration from the soil and epikarst zone, towards vadose zone and phreatic zone (Ghasemizadeh et al. 2012). Accordingly, all subsequent models adopted a similar divisional structure formed by four main components: the soil zone, epikarst zone, vadose/unsaturated and phreatic/saturated zones. Such models adequately describe the structure of karst systems formed by epigene processes.
In this paper, we focus primarily on karst systems formed by epigene and/or hypergene processes; we do not discuss hypogene karst whose genesis and development may be driven by very different mechanisms. In epigene karst systems, the hydrological interactions between meteoric water, surface water, soil water and karst groundwater are complex. During rainfall events, the meteoric water first enters the soil zone or/and the epikarst where soils are absent. Within this zone, rapid flow may occur in vertical fractures and conduits (sinkhole), while slower flow takes place in matrix. The epikarst zone can also play a major role in concentrating water flow and thus generating intense hydraulic response that may propagate through the system and be monitored at karst springs. In shallow karst system, the epikarst processes thus can dominate the spring discharge and chemical signal (Perrin et al. 2003). Although the intensity of fissures decreases rapidly with depth, concentrated water may continue to Fig. 1 Summary of karst modelling approaches that will be discussed in this review paper. Black font color indicates classical approaches; red font color indicates the approaches, variables or processes more recently considered in the modelling percolate quickly to the unsaturated zone through major tectonic fractures and shafts. Within this zone, a horizontal redistribution of the infiltrated water may become important. The unsaturated zone is mainly seen as a transmissive zone connecting the epikarst and phreatic zones, although it may be highly capacitive, especially when thick vadose zone exists and control the organization of the flows with the depth in the unsaturated zone of the karst system (Emblanch et al. 2003). Flow in the unsaturated zone can also dominate the transport processes and thus modify karst spring chemographs (Perrin et al. 2007). Once the water flows in the phreatic zone, the sub-horizontal flow component predominates and, although diffusive flow may still exist, groundwater is mainly transmitted by the network of karst conduits towards springs and outlets.
Another important characteristic of karst aquifer at larger time scale is that their hydrological functioning varies constantly due to the dissolution of carbonate rocks that continuously modify the properties of the flow path networks. The evolution process due to dissolution starts when meteoric water mixes with CO 2 in the soil zone. The majority of limestone dissolution occurs in the epikarst zone and to a lesser extent in the unsaturated zone and phreatic/saturated zone. Indeed, the water often reaches carbonate saturation when transmitted to the phreatic zone. The temporal variability of karst systems hydrological functioning is a key feature that distinguish them from classical aquifers where changes can occur over time, but not as much and as fast as in carbonate aquifers.

Modelling of key karst processes
The terminology of 'functional approach' and 'physically based approach' is considered when referring to the models of karst processes. In this section, we focus on the key karst processes that play important roles in the system's functioning: recharge, hydrodynamics, transport and speleogenesis. The commonly used methods for modelling these processes, either according to the 'functional approach' or 'physically based approach', are then reviewed.

Recharge process
Estimating groundwater recharge in karst is mostly relying on the modelling of the dynamics of the soil, epikarst and the unsaturated zone . The soil and epikarst have an important hydrological role, due to their high storage capacity and thus highly regulate aquifer recharge. After infiltrating through the epikarst, the recharge is generally divided into two major portions: the vadose percolation, which is the part of the diffusive flow through fissures and matrix, and the vadose fast flow, which is associated to piston movement of water flux through shafts, karstic conduits and fractures with large openings. Groundwater recharge can be modeled by means of diverse methods. The conventional karst recharge models are based on water balance or budget calculations (Blavoux et al. 1992;Scanlon et al. 2002). The methods rely on the fact that water entry should be equal to the amount discharged plus or minus the variation in the volume of water that is stored. By combining with geographical information system (GIS), an estimate of the spatial distribution of recharge can be obtained (Andreo et al. 2008;Malard et al. 2016;Radulovic et al. 2012). The water balance approaches have several advantages: they are easy to implement, but they suffer from several important drawbacks. First, they are based on empirical interpolations and estimates which may be erroneous. Also, they rely essentially on site-specific information, so that their transferability to other regions or towards the future is limited (Hartmann et al. 2014).
An alternative approach to karst recharge modelling relies on physically based distributed methods (Ireson and Butler 2011;Kiraly et al. 1995;Perrin et al. 2003). In this method, groundwater flow equations are solved under saturated or unsaturated conditions. Along this line, the pioneering work by Kiraly et al. (1995) considers a nested 2D model to simulate the epikarst, which redistribute flows toward a underlying 3D karst network model. Outflow from the epikarst is considered as recharge inflow at sinkholes or dolines and computed as a proportion of concentrated and diffuse recharge, respectively. The assumption of partial saturation is often adopted to model the flow in the saturated zone, which is under a steady-state regime. To deal with the condition of the total saturation, a time variant model may be used, but it requires much more data (Andreo et al. 2008). In this situation, recharge estimation with the functional approaches exhibit promising results and show good potential for estimating recharge and capture realistic flow behaviors in karst aquifers.
Karst recharge models can also be distinguished as a function of the considered approach, i.e. functional or physically based, accounting for unsaturated zone flow processes. In the case of limited data availability and/or aquifer scale modelling purposes, functional recharge models are generally preferred (Fleury et al. 2007;Jukić and Denić-Jukić 2009;Tritz et al. 2011). However, contrary to physically based approaches (Andreu et al. 2011;Dvory et al. 2016;Hughes et al. 2008;Ireson and Butler 2011;Kiraly et al. 1995;Perrin et al. 2003), such functional approaches do not allow considering the spatial variability of recharge and can only provide a mean recharge, spatially averaged over the karst catchment. In the case of physically based approaches, spatial variability of soil and epikarst hydrodynamic properties can be measured in the field and considered either for the explicit modelling of water movement in the vadose zone (Kordilla et al. 2012;Rimmer and Salingar 2006) or the estimation of groundwater recharge, distributed over the karst catchment.

Dual-medium flow
Karst aquifers functioning is often characterized by the dual flow behavior, which results from high-conductive flow paths made of karst conduits and fractures, within a low conductive matrix. Such conceptualization of flow is also well-suited to describe flow in either the epikarst or vadose zone. However, due to multiple phase flow processes (mainly air and water) in these compartments, it was mainly considered for the modelling of single-phase flow in the saturated zone. The first attempt to simulate groundwater flow and solute transport in a karst aquifer with a hybrid discretecontinuum numerical model was proposed by Kiraly (1998), who considered conduit flow in discrete pipes embedded in a continuum porous medium with Darcian flow. This hybrid approach was implemented for the development of the CAVE (Carbonate Aquifer Void Evolution) model (Clemens et al. 1996;Liedl et al. 2003).

Filled conduit flow/laminar and turbulent flows
Early flow models only accounted for laminar flow, although turbulent or free surface flow may occur in conduits or other enlarged flow paths. Turbulent processes were considered in MODFLOW-CFP (Hill et al. 2010;Reimann and Hill 2009) while considering the hydraulic gradient in conduits and adjusting the hydraulic conductivity so that it accounts for turbulence (Shoemaker et al. 2008) when the hydraulic gradient exceeds a critical value (Fig. 3a). Further developments were proposed by Reimann et al. (2011Reimann et al. ( , 2014 to simulate conduit-associated drainable storage and unsaturated conduit flow, and account for time variable boundary conditions (Fig. 3b). MODFLOW-2005 considers laminar flow in the continuum (calculated by Darcy's equation), while CFPM1 (Conduit Flow Process Mode 1) considers laminar flow in the continuum (calculated by Darcy's equation) but also laminar or turbulent flow in a discrete pipe network (calculated by the Hagen-Poiseuille and the Colebrook-White equations). CFPM2 (Conduit Flow Process Mode 2) considers both laminar flow (calculated by Darcy's equation) and turbulent flow in the continuum (estimated by adjusting the hydraulic conductivity as described by Shoemaker et al. 2008). Finally, CTFc (Conduit-Type Flow in Continuum Cells) also considers laminar flow (calculated by Darcy's equation) and turbulent flow in the continuum (estimated a power law).
The simulations proposed by Reimann et al. (2011) consider nonlinear flow at the outlet of a karst conduit embedded in a porous matrix, for a direct recharge pulse injected in the conduit and a constant diffuse recharge in the matrix. These simulations show that the non-consideration of turbulent flow processes may lead to an overestimation of the discharge at the outlet of the karst conduit (Fig. 3b).

Partially filled conduit flow/free-surface flow
The simulation of strongly variable transient flow processes, such as flood, requires considering free-surface flow, wave propagation in karst conduits, and changes between pressurized and non-pressurized conduit flow. The Saint-Venant equations for free-surface flow simulation in karst conduits were implemented in the coupled overland and groundwater flow model MODBRANCH, to simulate unsteady and nonuniform flow processes in karst conduits. The Saint-Venant equations allow describing unsteady flow in partially filled conduits (de Saint-Venant 1871) and are therefore often used to simulate discharge in conduit-flow-dominated karst systems. In the unsaturated zone, free-surface conditions control the conduit flow, so that Manning's equations are generally considered, although they may underestimate head losses due to turbulent flows (Jeannin and Maréchal 1995).
As the cross-sectional area of karst conduits is finite, both pressurized and non-pressurized flow may occur. Equations were therefore adapted for pressurized conditions through considering the Preissmann slot model in MODBRANCH package (Reimann et al. 2011), or the Darcy-Weisbach equation in pipe-flow models (Jeannin 2001), such as the Storm Water Management Model (SWMM) (Campbell and Sullivan 2002;Chen and Goldscheider 2014;Peterson and Wicks 2006). However, the use of such models requires detailed distributions of hydraulic parameters which are often unavailable due to the unknown location and geometry of karst conduits.

Transport processes
The highly heterogeneous flow fields in karst media originating from the varieties of flow paths can result in complex transport phenomena in the system (Fig. 4a, b). Indeed, the range of residence time of karst groundwater can be considerably contrasted, spanning from hours to decades. Moreover, the mean residence time typically exhibit a strong temporal variability and is highly dependent on the hydraulic condition, which greatly controls transport and speleogenesis processes. Accordingly, various computational models have been developed to simulate transport dynamics in karst formations (Fig. 4c, d).
Most interpretation methods of artificial tracer data are based on the classical one-dimensional advection-dispersion equation (ADE). Field observations have shown that the breakthrough curves from karst aquifer are most likely asymmetric and exhibit a sudden rising limb, a long tail, and multiple peaks Moreno and Tsang 1991). The traditional ADE performs poorly in characterizing these BTCs. To handle the tailing effect, several transport models have been proposed, which may generate more adjustable tails (Morales et al. 2010). The Mobile-Immobile Model (MIM) originally proposed by Coats and Smith (1964) partitions the aquifer into mobile and immobile regions. The method was later successfully applied to laboratory column experiment data by van Genuchten and Wierenga (1977) and field data (Atkinson et al. 1973;Brown and Ford 1971;Field and Pinsky 2000;Göppert and Goldscheider 2008;Massei et al. 2007Massei et al. , 2006Mull et al. 1988;Toride and Leij 1996), and incorporated into standard solute transport codes such as the popular MT3D (Zheng and Wang 1999). Skopp et al. (1981) proposed dual-permeability model which modifies the MIM by considering both regions are mobile and ADE is valid in both regions; this model is referred to as Multiple Region Advection-Dispersion (MRAD) by Majdalani et al. (2018). Another possibility for simulating the long tails of karst BTCs is the Continuous Time Random Walk (CTRW) method (Berkowitz et al. 2006). The CTRW approach visualizes solute transport as a series of particle jumps or transitions with spatially changing velocities. Recently, it has been applied to tracer test interpretation in the context of karst aquifers by Goeppert et al. (2020).
To handle the multi-peak effect, some progresses have also been made. The dual-peaked BTCs were successfully modeled by the application of the MIM (Field and Leij 2012; Maloszewski et al. 1994Maloszewski et al. , 1992Werner et al. 1997). Field and Leij (2012) compared the ability of dual-advection dispersion equation (DADE) and the weighted sum advection-dispersion equation (WSADE) to produce dual-peaked BTCs. Recently, Wang et al. (2020) applied the DADE model and successfully fitted the BTCs generated by controlled laboratory conduit network experiments.
Theoretical studies have been conducted to investigate the problem of solute transport in a system of a single conduit embedded in the background porous media (Li et al. 2008). Analytical solutions have been derived for these problems. Although the hydrodynamic interaction between the conduit and matrix is considered, the mass exchange at the conduit-matrix interface is assumed constant and does not vary according to the hydrodynamic condition (Li and Field 2016;Li and Loper 2011). Zoghbi and Basha (2020) extended this model to incorporate a fully coupled conduit-matrix exchange behavior, accounting for the change in hydraulic gradient/flow and transport between the conduit and matrix. The application of this model to field data obtained on two karst aquifers demonstrated its effectiveness in simulating karst BTCs and estimating key transport parameters.
When one attempts to model spatially variable output phenomena, physically based models may provide more accurate representation of some of the transport processes (Birk et al. 2006(Birk et al. , 2005Faulkner et al. 2009;Oehlmann et al. 2015). In these models, spatially variable velocities were estimated first from flow simulations and then used in the transport simulation. Thus, depending on the parameterization of flow models, the type of distributed transport model includes equivalent porous medium (Zheng and Wang 1999), dual-porosity (Gerke and Van Genuchten 1993), discrete conduit network (Peterson and Wicks 2006) and hybrid (Faulkner et al. 2009;Oehlmann et al. 2015) approaches.
Many characteristics/processes of karst transport dynamics have not been explicitly considered in the existing karst transport models. For instance, inertial flow and dynamic transition between laminar and turbulent flow regimes that frequently occur in karst conduits may have important effects on the transport process but are difficult to handle with traditional numerical models. To simulate these processes, computational fluid dynamics (CFD) simulations need to be performed (Hauns et al. 2001). However, the use of CFD models is limited for the understanding of general transport behavior based on scenario models, as they require significant computational resources. Therefore, their application for regional-scale karst problems is prohibitive.
The lattice Boltzmann methods (LBMs) offers an alternative numerical approach for simulating the effect of these complex hydrodynamics on conduit transport properties while consuming acceptable computational resources (Sukop et al. 2008;Anwar and Sukop 2009).

Speleogenesis
Earliest speleogenesis models were proposed by Dreybrodt (1990Dreybrodt ( , 1996 and Palmer (1991). These one-dimensional (1D) models allowed demonstrating the importance of nonlinear dissolution processes and the existence of positive feedback loops while considering the evolution of a single fracture (Fig. 5a). They also allowed calculating dissolution length scales and gave rise to the concept of breakthrough time; the breakthrough time corresponds to the time when ; t = 25τ, t = 50τ, t = 75τ, and t = 100τ (τ is the time it takes for the initial fracture aperture at the inlet to double); d concentration and flow discharge evolution during karst generation in a 20 m long fracture for linear and non-linear dissolution kinetics, in which q 0 is the initial flux in the fracture, lp is the penetration length, and C sat is the concentration at saturation (after Szymczak and Ladd 2011) a steep increase in the flow rate is observed, due to the increasing widening of the fracture aperture at the outlet Dreybrodt (1990). These models were the basis for more complex 2D models (Dreybrodt 2012; Groves and Howard 1994;Kaufmann et al 2012;Siemers and Dreybrodt 1998), with more complex dynamics related to 2D networks of fractures (Fig. 5b).In such cases, the region with dissolution of first-order kinetics penetrates downstream, as the discharge through the fracture network increases (Hubinger and Birk 2011), which results in a main branch of the karst conduits network fanning out at the discharge side (Fig. 5b).
At the scale of a single fracture, 2D models allowed demonstrating that aperture heterogeneity can be at the origin of preferential flow paths and thus generate earlier breakthrough than in the 1D case (Hanna and Rajaram 1998). Flow exchange between matrix and fracture (Bauer et al. 2003;Liedl et al. 2003;Kaufmann and Braun 2000), as well as flow exchange between prominent fractures and fissures (Gabrovsek et al. 2004), can also generate earlier breakthrough. Fingering of the dissolution front within a single fracture (Fig. 5c) was further demonstrated as the consequence of the high instability of the dissolution front propagation (Szymczak and Ladd 2011). In such modelling, breakthrough time and concentration profiles, as well as flow discharge evolution during speleogenesis highly depend on the considered dissolution kinetics (Fig. 5d). Further modelling was performed to study the influence of mixing corrosion and the competition between different flow paths and the effect of CO 2 sources (Gabrovšek et al. 2000), flank margin caves formation (Dreybrodt et al. 2009;Dreybrodt and Romanov 2007), karstification around dam sites (Dreybrodt et al. 2002), as well as buoyant convection (Chaudhuri et al. 2013). 2D fracture networks dissolution models allowed analyzing the horizontal and vertical evolution of caves (Gabrovšek and Dreybrodt 2001;Groves and Howard 1994;Siemers and Dreybrodt 1998). Dual-porosity models were further developed (Kaufmann and Braun 2000;Liedl et al. 2003), while considering flow through discrete conduits and through the porous rock matrix as well as exchange between both media. 3D model accounting for both speleogenesis and landscape evolution was proposed by Kaufmann (2009) to assess karst evolution over time. A complete review of fracture network speleogenesis models has been proposed by Dreybrodt et al. (2005) and Dreybrodt and Gabrovšek (2018).
The topology of the initial flow path network, which often consists of discrete fracture networks, may also influence significantly the dissolution pattern (Aliouache et al. 2019). When the flow is parallel to the persisting fracture orientation and becomes channelized (Fig. 6a-c), the initial hydrogeological setting highly controls the dissolution pattern. On the contrary, intermittent interruptions of the positive feedback loop between flow and dissolution occur when the flow direction is normal to the persisting fracture orientation (Fig. 6d-f), and the breakthrough is delayed.

Summary with a view towards hydrogeological modelling
From the above discussion, it can be seen that the modelling of key karst processes can be achieved by either functional or physically based approaches, which have important differences in the conceptualization of karst systems. Table 1 presents a comparison of the functional and physically based models. The functional modelling approach allows bypassing the geometrical complexity using equivalent behavioral parameters that reflect the overall hydrologic and hydrodynamic functioning of the karst system. On the other hand, physically based modelling scheme treats the system as an assemblage of interacting individual components and allows for the explicit integration of complex physical laws for hydraulic processes and interactions. However, some of the model parameters (e.g., conduit friction coefficient, exchange coefficient) may need to be obtained by numerical calibration rather than from direct measurements from the field. Furthermore, the computational time for solving physically based models can be substantially larger than that for functional models.
Since the important characteristics of hydraulic processes are preserved in the conventional functional and physically based modelling, some hydrogeological characteristics of karst aquifers can be reproduced. However, simplifications are commonly adopted in order for the problem to be tractable. Many aspects may be improved for a better modelling of karst hydrological functioning. To improve recharge simulation methods, a promising direction for future research is the development of hybrid models that assimilate the advantages of both functional and physically based models. Some of the current recharge estimation models have already exhibited some of such characteristics. For example, Bradley et al. (2010) combined several flow paths and storages within a functional model. In the subsequent stochastic recharge model proposed by Hartmann et al. (2012Hartmann et al. ( , 2015, statistical distribution functions were used to represent the soil and epikarst heterogeneities controlling storage and flow paths. Most physically based models that were developed in the past focused on the behavior of the karst aquifer in the phreatic zone. Despite recent progress in the identification of karst vadose zone dynamics, the appropriate characterization of process parameters and related conceptual models remain to be initiated for an efficient implementation of the typical processes of karst vadose zone into respective physically based models. Complete models should therefore account for unsaturated flow in the vadose zone, which would require considering film flow in wide fractures and shafts, and saturated storm flow in smaller fractures. These models should also consider the dampening effect of the epikarst related to temporary storage and retardation of precipitation. In general, the simulation of coupled saturated and unsaturated flow is still a challenge in hydrogeology, especially for karst systems (Reimann et al. 2011;de Rooij et al. 2013). The dual continuum approach may be adopted for such a purpose Teutsch and Sauter 1991). However, a sacrifice for reality of aquifer geometry must be accepted (Kordilla et al. 2012). A promising way of addressing the difficulties in modelling complex karst hydraulics is the employment of Lattice Boltzmann Methods (LBMs). The LBM-based karst models are particularly attractive as they are capable of simulating inertial flows, incorporating complex conduit geometries, and solving coupled Darcian and  Navier-Stokes flow problems. They also allow solving solute transport problems characterized by laminar and eddy flows in conduits or Darcian flow in porous matrix, and tackling multiphase flow and transport problems (Anwar and Sukop 2009;Sukop et al. 2008).

Modelling of karst conduits network
In karst aquifers, the majority of flow is conducted through the conduit network. The distribution of conduits within the aquifer influences the hydrodynamic functioning of karst aquifer. Thus, they should be ideally incorporated into numerical models, particularly when the physical mechanism of flow and transport in karst systems is concerned. The rigorous analysis of parameter uncertainty requires a large number of model-runs with an ensemble of plausible karst conduit networks. Therefore, methods for generating realistic karst conduit networks are needed. In practice, karst networks may be generated (1) based on field-mapped information of conduit geometry; (2) by stochastic simulation or (3) by sophisticated transport models accounting for complex geochemical reactions and dissolution processes. In the literature of karst hydrology, the stochastic and geochemical approaches are also referred to as structure-imitating and process-imitating approaches, following the classification of numerical methods proposed originally by Koltermann and Gorelick (1996) for simulating heterogeneity in sedimentary deposits (Borghi et al. 2012;Pardo-Iguzquiza et al. 2011). The structure-imitating approaches relies on a statistical approach aiming to reproduce the structure of conduit network without considering physical and chemical processes, while the process-imitating approaches rely on dissolution processes that control conduit network evolution, i.e., speleogenesis.

Field-mapped (and interpolated) karst conduits networks
Conduit profile and pattern may be determined from detailed speleological investigations (e.g., Jouves et al. 2017) which allow reconstructing the three-dimensional geometry of explored cave systems. During such a field survey, threedimensional information may be acquired at a sequence of topographic stations before being processed numerically to obtain discrete conduits with elliptical or rectangular section shapes. A 3D conduit network may then be constructed/ simulated by interpolating the two-dimensional sections. This task can be solved with common speleological programs such as Survex, Visual Topo or GHTopo. Rendering with more realistic section shapes is also available in these programs. LiDAR technologies also allow to precisely map karst conduits and caves (Jaillet et al. 2011). This method is adapted to explore individual conduits.
As only a small portion of conduits can be surveyed, modelling techniques are required to reconstruct 3D karst conduit networks based on one-and two-dimensional measurements. Henrion et al. (2010) proposed a new method (object-distance simulation method-ODSIM) to generate cylindrical and lens shaped envelops along an irregular curve considered as the skeleton of the karst conduit, in order to reproduce its 3D geometry. Rongier et al. (2014) further improved this method by considering a fast-marching method to account for the different geological properties influencing conduit shapes (Fig. 7).
The method produced conduit network with more realistic section shapes. The use of these generated complex 3D conduit network for flow and transport simulation at aquifer scale remains to be demonstrated, considering the high computational cost of 3D simulations to couple geometry modelling with modern simulators of flow and transport.

Stochastically-generated karst conduits networks
The complete geometrical, geomorphological and topological properties of karst systems cannot be fully assessed due to the difficulty to identify karst conduit networks at all scales. This favored the development of stochastic approaches that rely on the hierarchy and scaling of natural karst systems, which can be adequately characterized by fractal geometry or power law models (Hendrick and Renard 2016a;Pardo-Iguzquiza et al. 2011).
The spatial organization of natural karst networks, which often exhibit a hierarchical structure as highlighted by extensive field observations, can be characterized by a fractal dimension, D, which quantifies the manner of how fractals spread and cluster in the Euclidean space. The reported fractal dimension for natural karst networks range from 1.043 to 1.67 (Hendrick and Renard 2016b;Pardo-Igúzquiza et al. 2012).
The fractal characteristics may be incorporated in stochastically-generated karst networks. One example is the diffusion-limited aggregation approach of Pardo-Iguzquiza et al. (2012). In this approach, the geometry of conduit segments (e.g., orientation and length) is first simulated by sampling with replacement from known conduits (conduit templates). Then, the topology (connectivity between conduit segments) is generated stochastically using the diffusionlimited aggregation method. Although the method is able to incorporate the fractal character of karst conduit geometry, the generated connectivity is not explicitly based on statistical input parameters. Thus, it is not guaranteed that the resulting network will have a realistic connectivity pattern, even though this approach was successfully applied for either conditional or non-conditional stochastic simulation of karst conduit networks.
Apart from the fractal approach, other numerical methods may be used to simulate karst. Borghi et al. (2012) used a fast-marching algorithm to simulate the evolution of complex karst patterns formed under different boundary conditions (Fig. 8a-c). Jaquet et al. (2004) developed a latticegas cellular automation approach to simulate the diameter growth of a conduit network and the spatiotemporal evolution of a multiphase speleogenesis system (Fig. 8d). Ronayne (2013) used a non-looping invasion percolation model (Stark 1991) that allows generating dendritic (i.e., non-looping) connectivity patterns. Although dendritic geometries are generally related to stream networks, such connectivity patterns can be generated by a focused recharge, which results in dendritic karst conduit networks (Palmer 1991).

Geochemically-simulated karst conduits networks
Geochemically based models, which benefited from the increasing knowledge of speleogenesis processes, consider carbonate dissolution kinetics to simulate conduit growth and karst network evolution in response to hydrologic and geomorphic variations. Such model can be used to reproduce natural karst patterns in response to paleo-hydrological condition. One of the approaches is the hydrochemical model of de Rooij and Graham (2017), which has been applied by Henson et al. (2018) to understand the major factors, which contribute to the formation of large springs. However, such methods favor the development of hierarchical networks, which are one particular geometry among many other possible geometries of karst networks (Palmer 1991). The computational cost of these speleogenesis models is also very high and they are very difficult to constrain with field data.

Summary with a view towards hydrogeological modelling
Each of the three types of conduit network generation models have their strengths and limitations (Table 2). Fieldmapped karst networks models allow capturing many realistic features of karst conduits; however these models become limited to characterize phreatic conduit network in deep rocks and caves with small diameter that human cannot access. The stochastic conduit network approaches have the merit to be efficient and also applicable for 3D problems. However, they require a detailed statistical analysis of conduit properties and accurately deriving control parameters (e.g., conduit sizes and fractal dimension). Such a practice requires a large amount of data at all scales, and is therefore barely difficult to apply for realistic problems. The processimitating karst models, which may capture some genesis characteristics of natural karst, are sensitive to the inferred paleo-flow fields and boundary conditions. These models need to be further improved to take into account geological, thermal and geomechanical effects to reproduce actual karst systems. Recent advances in karst modelling have proposed such evolution but still remain unapplied to realistic karst system.
To generate better karstic conduit networks, a future research direction that is attracting much effort is the development of process-oriented numerical models for catchment scale problems. Several studies adopted such a modelling strategy and have already shown promising results. For example, de Rooij and Graham (2017) used the spring locations to condition a 2D conduit network generated by a hydrochemical model. Aliouache et al. (2019) applied numerical karst genesis modelling in a geologically-mapped fracture network following dissolution kinetic law. Their conduit growth is regulated by the realistic topological properties of natural fractures. Jiang et al. (2022) proposed a hydro-chemical model to simulate incipient karst generation in 3D jointed layered carbonates and highlighted the role of different flow regimes as well as the initial penetration length in the generation of distinct karst morphologies (i.e. pipe-like, stripe-like and sheet-like karst features). Such karst genesis models may thus supply additional information for outlining the geometry of conduit network in realistic modelling practice.

Modelling of karst hydrogeological functioning
Various analytical and numerical approaches have been used to model the hydrogeological functioning of karst aquifers. The preference for a functional or physically based modelling scheme often depends on the objective of modelling application, scale of problem, system complexity and data availability. The functional approach  Borghi et al. 2012)) and with the cellular automation approach (d after Jaquet et al. 2004) has the advantage of simplicity and its great efficiency to model large-scale system functioning. In contrast, the physically based models can straightforwardly incorporate spatial heterogeneities and model complex hydrodynamic processes. Depending on the discretization type, the physically based models may be further classified as equivalent porous medium (EPM), double continuum (DC), discrete conduit network (DCN) and combined discrete-continuum (CDC) models as proposed in former reviews (Hartmann et al. 2014;Kovács and Sauter 2007). In this section, commonly used models for simulating hydrogeological behaviors of karst systems, including spring responses to natural forcing (i.e., rainfall or recharge) and borehole responses to artificial forcing (i.e., hydraulic tests) will be reviewed.

Karst spring responses
One peculiarity of karst aquifers is that most of the groundwater discharge occurs through a small number of outlets, and sometimes through a single spring, which can be characterized by hydrographs or chemographs controlled by the whole water moving through the aquifer. Numerical studies therefore focused on the simulation of karst spring responses, with either functional or physically based model, to understand system's dynamics.

Functional approach
Functional approaches model the karst system's functioning or evolutionary dynamics at the scale of the entire system, from the infiltration of meteoric water into the soil zone and epikarst to discharge at springs. The dynamic processes, such as infiltration and storage in various aquifer components are modeled with linear or nonlinear empirical relationships. Due to this lumped nature of model structure, the parameters of the functional models only represent effective properties of the whole karst system and may not bear a physical meaning. For this reason, the model parameters cannot be quantified by field measurements and only the dynamic of discharge at springs can give insight about the base flow recession that can be further considered in the model. Instead, such parameters may need to be determined by an indirect process of numerical calibration. Calibrations of functional karst models are important task, particularly when they are used for real applications. The commonly used data for calibrating functional karst models include continuously monitored flow and physicochemical data at karst springs and to a much less extent, in-situ measurements of water levels in monitoring boreholes (Ladouche et al. 2014;Mazzilli et al. 2019Mazzilli et al. , 2011. Functional models can be applied for single event (or recession) analysis or time series analysis. Recession analysis is widely used to understand the hydrodynamic functioning of karst system (e.g., Amit et al. 2002;Bonacci 1993;Chang et al. 2015;Cinkus et al. 2021;Fiorillo 2011;Florea and Vacher 2006;Forkasiewicz and Paloc 1967;Kovács et al. 2005;Malík and Vojtková, 2012;Mangin 1975 Rahnemaei et al. 2005). Detailed mathematical analysis has been continuously and progressively inserted to increase the efficiency of this method (Gárfias-Soliz et al. 2010;Jukić and Denić-Jukić, 2011;Mathevet et al. 2004). The early stage of karst functional models focused on reproducing the discharge behavior caused by the duality of flow regime, one of the most important concepts in karst hydrology. The components of fast conduit flow and diffusive matrix flow are simulated using two separate reservoirs or tanks, which are combined to generate the discharge. These models were extended to consider the conduit-matrix coupling by formulating an exchange relationship between the reservoirs. For instance, in the models proposed by Geyer et al. (2008), matrix may recharge the conduit reservoir if a pre-defined hydraulic threshold is reached. Functional models that include 'two-way' flow exchanges between matrix and conduits have also been proposed (Mazzilli et al. 2011). Another alternative way of simulating the dual flow characteristics of karst aquifer is to couple a matrix reservoir and a routing function representing the rapid flow component, which is transferred directly to spring discharge (Le Moine et al. 2008).
To account for the influence of soil and epikarst in regulating recharge and providing storage, a separation function that divides recharge into matrix and conduit reservoirs according to a pre-defined fraction ratio has been introduced in classical functional models (Charlier et al. 2012;Fleury et al. 2009). Another approach developed to account for the particular flow dynamics related to soil and epikarst consists in integrating a hysteresis scheme that allows activating and deactivating concentrated recharge (Tritz et al. 2011).
A new class of functional karst models that integrate the hydrochemical modelling into original hydrodynamic analyses have also emerged. Thorough incorporating extra hydrochemical data, such as chemical and heat tracers (Anderson 2005;Covington et al. 2011;Luhmann et al. 2011;Martin and Dean 1999;Saar 2011), water chemistry (Barberá and Andreo 2012;Cholet et al. 2015;Filippini et al. 2018;Mudarra et al. 2014;Raeisi and Karami 2011;Toran et al. 2006), water conductivity (Massei et al. 2006;Wang et al. 2020) in addition to hydrographs, it has been shown that these integrated functional models provides a more reliable characterization of the system functioning of karst aquifers. Hydrochemical information in functional models (Charlier et al. 2012;Hartmann 2013;Hartmann et al. 2016Hartmann et al. , 2017Sivelle et al. 2022) were introduced to combine information inferred from both hydrochemical signal and hydrodynamic signal.
The simplicity of lumped model structure may be regarded as an advantageous aspect, because detailed information, such as geometrical and hydraulic properties for dominant processes, is often not available when analyzing complex karst systems, for which physically based modelling is greatly limited. However, it is still very important to continue improving the realism and accuracy of functional karst models, because the spring responses (i.e., the spatial and temporal variations of discharge and physicochemical parameters) are governed by the recharge pattern and physical processes integrated within the model. Developments are needed towards (1) a more thorough characterization of the flow and transport dynamics respecting more details of real karst system, such as the complex role of epikarst storage (Perrin et al. 2003) in the regulation of recharge, mixing of waters from rainfall and aquifer tributaries, as well as chemical reaction, and (2) a more rational and efficient implementation of model structure respecting the identified dominant processes and accounting for the coupling between some of these processes. For more extensive discussion on the use of functional models for modelling karst spring responses and other improved implementations, the reader is referred to Jeannin and Sauter (1998), Kovács and Sauter (2007) and Hartmann et al.(2014).

Physically based approach
Extensive studies have been conducted to model karst spring dynamics and understand their relation to spatial heterogeneities and physical processes (Doummar et al. 2012;Kordilla et al. 2012;Liedl et al. 2003;Oehlmann et al. 2015;Xanke et al. 2016). The increased knowledge of the mechanisms of conduit genesis and hydraulics promoted the development of physically based karst models that incorporate the mechanisms governing the flow and transport processes and simulate karst spring signal evolution as a hydrological response to spatial heterogeneity.
Various numerical discretization methods have been proposed and the one based on equivalent porous medium (EPM) concept is frequently used (Ghasemizadeh et al. 2015). In the EPM model, the rock matrix, fissures and conduit networks are represented by a continuum domain with equivalent hydraulic conductivities assigned to individual grid cells. The inclusion of conduit features can be achieved through using lines of drain cells and internal sinks (Quinn et al. 2006;Quinn and Tomasko 2000) or lines of high hydraulic conductivity cells (Lindgren et al. 2005). Since the modelling objective time dependent behaviors such as spring discharge and water table variation, transient simulations are generally required. Storage parameters, such as specific yield, may be assessed by core testing techniques or in some studies were treated as variables to be determined by a calibration process. The transient recharge is difficult to quantify as it depends on various meteoric factors such as rainfall and evapotranspiration as well as reservoir characteristics, such as elevation, runoff and recharge area. Water budget methods are often adopted to identify the recharge type and estimate the recharge parameter values. The initial and boundary conditions of the karst groundwater system may be identified and approximated by hydrological and hydrogeological characterization of field observations and monitored data (Ghasemizadeh et al. 2012).
In general, the EPM models are more suitable for karst modelling at regional scales where the salient influence of conduit hydraulics is averaged out. Scanlon et al. (2003) developed an EPM model for simulating regional scale groundwater flow in the Edwards aquifer in Austin, Texas and found that the model was able to simulate temporal variations in spring discharge and to reproduce the effect of pumping on discharge with fairly good accuracy. To evaluate the potential impact of conduit flow, Saller et al. (2013) integrated high conductivity zones, which represents karst conduits, into the EPM models developed originally by Long and Putnam (2009) for the Madison aquifer. They observed that the inclusion of high conductivity zones resulted in a better modelling of transient water table but did not improve the capacity of the original EPM model in karst spring discharge simulations. The classic physically based numerical models based on the EPM formulation can capture some aspects of the karst hydrological behavior, but were inadequate to simulate solute transport. Worthington and Ford (2009) created a regional model for the limestone aquifer at Mammoth Cave (Kentucky, USA) and investigated the difference in model response with and without the presence of high-permeability conduits. The authors found that the EPM model, even while considering conduit effects, cannot reproduce the tracer test results.
Apart from the EPM approach, discrete conduit network (DCN) model was used to successfully simulate flow of a karst system in Kentucky, USA (Thrailkill 1974). Jeannin (2001) modelled the hydraulic response of the Holloch cave (Switzerland) with a DCN model, and successfully related it to the local cave geometry. In subsequent studies, the storm water management model (SWMM) that is based on the DCN formulation has been adopted for simulating flow and transport in conduit-dominated karst systems (Campbell and Sullivan 1999;Chen and Goldscheider 2014;Peterson and Wicks 2006;Wu et al. 2008). One of the main focuses of such models was the examination of the influence of conduit hydraulic properties. Peterson and Wicks (2006) performed a parameter sensitivity test using SWMM, and reported that a slight change in the relative conduit roughness resulted in a substantial variation in the simulated spring discharge and solute transport. Because SWMM solve one-dimensional flow and transport equations, the calculation is very efficient even for solving regional problem. It can be coupled to a calibration system to study model parameter sensitivity and fitted to field measurements of long-term time series such as spring discharge and hydrochemical variation (Chen and Goldscheider 2014). Such a calibration capability also offers a possibility of inferring essential geometrical structure of karst network through inversions of multi-tracer test including transient flow (Sivelle et al. 2020). The DCN methods focus on representing conduit hydraulics but do not address the flow and transport processes in the matrix adequately.
To overcome such a deficiency, a dual continuum model may be adopted. The dual continuum approach is able to simulate the dual flow behavior of karst system without the need to incorporate an explicit representation of complex conduit networks (Teutsch and Sauter 1991). The integration of saturated and unsaturated flow regimes with the dual continuum conceptualization has also been implemented in the HydroGeoSphere model (Therrien et al. 2010).

Hybrid approach
The strong hydraulic coupling effect of various zones in a karst system can affect significantly the simulated aquifer functioning behavior. For example, the recharge function has a strong impact on spring hydrographs in karst aquifers (Hartmann 2013;Kiraly et al. 1995;Kovács et al. 2005). One promising method is thus to include recharge through the use of a functional reservoir model. In this way, the role of the epikarst zone in regulating and storing rainfalls is considered and the output of the reservoir model gives the recharge input for the physically based model component (Borghi et al. 2016;Chen and Goldscheider 2014). This method can be implemented in all kinds of karst simulation model and has a potential benefit to simulate a real karst aquifer. Ofterdinger et al. (2015) extended the simple hybrid models and used a complex hydrological model to compute recharge distribution and the runoff of precipitation. However, this kind of model is very difficult to parameterize given that the number of measurement stations is often limited. On the other hand, the parameter values of recharge models are difficult to identify. One encouraging way to address these difficulties is to incorporate some of the recharge model parameters into an inverse framework and to estimate their values by an optimization process (Mohammadi et al. 2018;Kavousi et al. 2020). The calibrated parameters need to be back-checked with field data to confirm if they are physically meaningful. Thus, their application to real problems is only restricted to well monitored karst systems.

Numerical representation of hydraulic tests
Karst aquifer characterization is typically based on spring response interpretation. Recently, it has been extended by systematic analysis of pumping tests. Experimental field sites equipped with boreholes and in-situ monitoring sensors may provide useful information about karst aquifer properties at local scales. However, they do not allow determining hydrodynamic properties at catchment scale because of the restricted extent of the investigated area due to limited pumping rates. In particular cases with exceptional pumping rate of several hundred liters per second, large-scale characterization of karst aquifers has been made possible (Maréchal et al. 2008;Noushabadi et al. 2011;Dausse et al. 2019).
Numerical models with distributed hydrodynamic parameter fields in both matrix and karst conduits can also be used to characterize karst aquifers at large scale. Discrete-continuum models represent a suitable approach to simulate karst aquifers (e.g., Kaufmann 2009;Kiraly 1998Kiraly , 2003Sauter et al. 2006) while considering the matrix interaction and highly anisotropic hydraulic parameter fields. Liedl et al. (2003) proposed a model comprising both a discrete-continuum for the simulation of laminar and turbulent pipe flow, and a continuum that represents the matrix. As this model applies steady-state pipe flow equations, it cannot capture the transient water storage within the conduit. To overcome this deficiency in representing realistic storage dynamics, the Saint-Venant equations for variably filled conduit flow can be solved (Reimann et al. 2014). Other important extensions for better reproducing the realistic hydrodynamic process in karstic media include coupling surface flow (de Rooij 2008) and variably saturated matrix flow (Rooij et al. 2013). Giese et al. (2017) used the discrete conduit-continuum model of Reimann et al. (2014) to investigate the effect of different boundary conditions as well as sink/source terms for idealized two-dimensional mixed karst aquifer systems. Giese et al. (2018) applied the same model and examined the effect of conduit flow regimes (i.e., laminar and non-laminar flow) on the drawdown in conduit and matrix with varied conduit hydraulic properties (i.e., conduit diameter and relative roughness). They found that the turbulent flow may generate important influence on pumping test interpretation for small conduit diameters, while for mature conduit system where the diameter is sufficiently large, the laminar flow equation may be applicable for simulating karstic conduit hydrodynamics even under the condition of fully-developed turbulent flow.
The numerical representation of pumping tests together with the borehole hydraulic head measurements can be used to infer the spatial distribution and hydraulic properties of karst conduit networks. Such an inversion approach is referred to as hydraulic tomography, which can be used to identify hydraulic connectivity of karst aquifers. While considering a quasi-Newton inverse method, Wang et al. (2016) estimated the transmissivity field and thus the major karst conduits and their connectivity at the scale of an experimental field site, which however highly depends on a priori information provided on the inversions (Fig. 9a, b). A stochastic Newton inverse method was further applied by these authors to quantify the uncertainty in the inverted transmissivity fields and address the issue of uniqueness of the inversion (Wang et al. 2017). On the same experimental field site, the cellular automata-based deterministic inversion Fig. 9 Karst conduit networks constructed by simultaneously inverting multiple cross-hole pumping test data; a inversion results with a constant initial hydraulic conductivity distribution; b inversion results with an initial hydraulic conductivity distribution inferred from interference tests; c inversion results obtained by the CADI method (after Fischer et al. 2017;Wang et al. 2016) (CADI) method (Fischer et al. 2017) allowed identifying a high transmissivity discrete conduit network embedded within a lower transmissivity matrix (Fig. 9c). The CADI method minimizes the difference between the observed and modeled hydraulic data by distributing the hydraulic properties along linear structures and iteratively modifying the conduit network geometry.
Further research considered harmonic pumping technique to assess both cross-borehole connectivity and karst conduit networks geometry (Fischer et al. 2018). With this technique, one can assess the degree of connectivity between the pumping well and the observation wells, which depends on the phase offset of the monitored responses, and the conductivity of main flow paths, which is controlled by the amplitude of the monitored response. High conductivity flow paths between observation wells and the pumping well can be identified with high frequency pumping, while dual connectivity between the observation wells and the pumping well (propagation occurring both in the karst conduit networks and in the matrix) can be highlighted by low-frequency pumping. The hierarchical flow path network (i.e., karst conduit, fracture, and matrix) can thus be assessed while considering various harmonic pumping tests frequencies. Tomographic inversions of harmonic pumping tests are also much faster (about 10 times) than the one using constant pumping rates and they require at least two times as few tests, which makes this method a powerful approach to enhance the interpretation of karst systems.

Comments on the modelling of karst hydrogeological functioning
In the conventional hydrogeological modelling, the karst is regarded as a particular and distinct type of aquifer. Besides, most of these hydrogeological modelling studies focused on the phreatic zone, with flow and transport dynamics in the epikarst and vadose zone either completely ignored or greatly simplified. Such treatments may be acceptable for modelling karst systems with thin and relatively homogeneous vadose zone, for which the effect of flow redistribution and storage of the epikarst and vadose zone may be negligible or approximated by simple recharge function. However, for most karst systems, the vadose zone is quite thick and affects the storage and recharge processes greatly. Under high flow conditions, fluid transfers are likely to occur at the conduit-matrix interface driven by high local pressure gradients. The matrix of the epikarst and vadose zone may thus be filled by conduit flows; while under low flow conditions, water stored in the epikarst and vadose matrix may be released to conduits and thus become the primary source of spring discharge. This regulation process can create retention and releasing of fluid and contaminant and thus affect greatly the hydrological and hydrochemical properties of karst aquifers and groundwater, as observed in the field (Perrin et al. 2003). Furthermore, the irregular roughness of karst conduits, which are often assumed to be straight line (in 2D) or cylindrical pipe (in 3D), may accommodate enhanced pressure loss under high flow rates, particularly at sites of diameter contractions, which has important effects on flow and transport processes. Hence, towards an improved hydrogeological modelling of karst aquifers, the capability of precisely modelling the coupled saturated-unsaturated hydraulic behavior is prerequisite, for which the important characteristics of flow and transport processes, e.g., conduitmatrix interaction, flow regime transition, need to be appropriately treated.
The development of calibration and validation procedures is a crucial issue for physically baased karst models based on field measurements, especially when they are applied for practical purposes. The genetically-generated conduit networks can be examined by comparing the geometrical properties of a simulated conduit network with those of the conduit pattern mapped by speleology exploration. Boreholes and underground research laboratories (Emblanch et al. 2003;Perineau et al. 2011) provide also valuable information on karst conduit sizes that intersect borehole and tunnel walls. These data might be of major importance for the calibration of karst network geometries. Hydrodynamic information, such as water table variations and pressure derivatives in response to pumping (single-well and cross-wells tests), as well as chemical concentration monitored during tracer tests (Ronayne 2013) can also help tuning the hydrological performance of karst conduit network models. Furthermore, the continuous physicochemical data monitored at springs and water table variations recorded in boreholes offer useful references for the calibration of karst models used for hydrological simulation at catchment scale. The newly-developed hydraulic tomography technique can be employed to calculate best-fit hydraulic diameter distributions for karst flow models based on an optimization process using multiple objective functions.

Outstanding issues and open questions
Modelling the hydrogeological and geochemical characteristics of karst aquifers is a challenging issue. In previous modelling studies, a reasonable understanding of the main components of karst aquifers and their hydrogeological functions has been established and efforts in developing quantitative techniques for modelling key dynamic processes have shown promising results. However, many important issues and scientific questions still exist. These issues/scientific questions, that would be valuable to investigate, offer opportunities for future research in the field of karst modelling.

Issues of functional karst models
The functional models based on empirical relationships have shown great adaptability and strength in simulating the hydrogeological behavior under different contexts and scales (Hartmann et al. 2015). However, it remains relevant to continue improving their realism and accuracy. Enhancements to the current modelling schemes by assimilating distinctive karst processes are worthy of further exploration. Some functional models have been proposed to condition the model to flow measurements in conduit or matrix (Fleury et al. 2009;Schmidt et al. 2014), to discharge at multiple springs (Rimmer and Salingar 2006), but also to consider the variations of recharge and storage in soil and epikarst in time (Tritz et al. 2011) and space . Understanding the sensitivity and better assess the parametric uncertainty of models (Cousquer and Jourde 2022) or their structural uncertainty is also an important aspect for future research. The predictive capacities of the functional models can be increased if field measurements in addition to spring discharge data are incorporated into the calibration procedure. Some studies have already been attempted to test the feasibility and consequence of integrating hydrochemical (Charlier et al. 2012;Hartmann 2013;Hartmann et al. 2016;Pinault et al. 2001) and gravimetric (Mazzilli et al. 2013) information. The joint use of data from different type of monitoring should help further reducing the ambiguity of functional models. Future studies should go on this direction.

Issues of physically baased karst models
In previous studies, simulation results obtained from physically baased karst models have shown consistency with well-controlled laboratory experiments and/or field tests. However, there remains a number of difficulties in applying physically baased numerical approaches for modelling karst processes. And a number of important issues still exist in numerical modelling of the hydrogeological behavior of karst aquifers using physically based models.
The first task will be the development of realistic hydrodynamic models that can represent the important transient flow and transport characteristics of natural karst systems (e.g., coupled free surface-Darcian flow systems, transition between saturated and unsaturated flow regimes, dynamic conduit-matrix interactions), and that can be used afterwards to capture the key hydrogeological and hydrochemical characteristics of karst systems (e.g., dual flow regimes, hydraulic thresholds, matrix retention and releasing of solutes).
The geometry of the created karst conduit networks needs to be consistent with the speleogenesis principles and field mapped conduit statistics so that the hierarchically-organized fluid flow can be properly modelled. Developments are needed towards (1) a more detailed characterization of the underlying statistics, such as multifractals for which more than one scaling exponent is expected (Hendrick and Renard 2016b;Pardo-Iguzquiza et al. 2011); (2) a more precise and efficient model to generate karst conduit network with more details to reproduce the diversity of individual conduit shapes and morphology (Rongier et al. 2014), the topological complexity in karst conduit populations (Pardo-Igúzquiza et al. 2019) and the correlation between conduit geometrical attributes (Jouves et al. 2017); (3) extrapolation of current 2D karst networks generated by the stochastic methods to 3D while preserving geostatistical, topological and hydrological properties.
Another major issue is to propose adequate upscaling approaches to hydrodynamic models to evaluate the largescale behavior of karst systems. The challenge for such models will be to preserve the hydrological and hydrochemical characteristics observed at the local scale. Some techniques were used to construct large-scale models on the basis of either upscaled permeabilities of local cells (Popov et al. 2009) or multiscale integration methods (Correia et al. 2020). In the future, more efforts should be put in other aspects, such as the development of more advanced parallelization schemes (De Rooij 2019), the simulation of multiphase flow effects on transport processes and, importantly, extension of the models to complex 3D systems.

General open questions: linking functional and physically baased models
As can be seen from above, the hydrogeological modelling of karst aquifers can be achieved by functional or physically baased approaches, which differ significantly in the conceptualization of geological media. Because the functional models treat the system as a black box and attempts to simulate the hydrogeological processes using analytical relationships between input signal and output response, the computational demand is much lower comparing to physically baased models. They may be easily coupled to an optimization system to calibrate model parameters and predict future aquifer's functioning behavior. However, an obvious drawback is that the calibrated model parameters may not bear any physical meaning. The continuum modelling scheme (i.e., EPM) mainly reflects the dynamic process in a karst from a more overarching view and simplifies the geometrical complexity using equivalent material properties derived from upscaling techniques. Such a treatment allows the continuum models to deal with large-scale problems with great efficiency. However, it cannot consider the effects of conduit flow dynamics, sharp hydraulic contrast in conduit and matrix, complex conduit-matrix interaction, and conduit surface dissolution kinetics.
On the other hand, the physically baased models regard the system as a group of interacting individual flow components and permits the integration of complex spatial variability for rock materials and realistic physical processes. Various physically baased models have been proposed to handle complexities at different level. A physically baased model can be established at a specific investigation scale without presuming the form of recharge-discharge transfer relationships. However, the resolution of discontinuous problems may require considerably larger computational time than the one needed for continuum models, which greatly limits the use of such models in calibration and inversion studies.
A direction for future research is to develop an integrated modelling approach by further exploring the link between functional and physically baased models. In practical applications, physically baased models that integrate all available data about the modelling site could be established first. It may then be used to drive functional proxies for the physically baased models, where the actual flow and transport physics is replaced by a set of reservoir structures characterized by empirical relationships that produce similar system responses. As such, the numerical efficiency will be greatly improved and the functional proxies could be adopted for calibration and sampling purposes at a much lower computational cost. To explore this link, a comparison of a large number of simulations needs to be conducted. For the task to be successful, advanced computing techniques, such as parallelization on GPU and/or CPU, as well as modern data analysis schemes, such as machine-learning and deep learning algorithms, may need to be implemented.
The improved understanding of the link between functional and physically baased model parameters may also help to propose more robust karst typology. Efforts to classification of karst system have been attempted previously with only functional approaches. The main focus was on looking for similarities in discharge responses. Although distinct discharge signatures have been recognized, difficulties still exist in extrapolating the results because the system's signature cannot be related to system properties (i.e., geometrical and hydrodynamic properties) which are not considered in functional models; as a consequence, the karst topology is not generalized. The use of physically baased models to infer structure and parameters of functional models will infuse physical meanings into the individual or groups of parameters. The development of such a physically baased karst typology is important for decision makings at the global scale under constraints of climate change and anthropogenic forcing as well as to improve our general understanding of karst hydrology.

Concluding remarks
This review paper began by presenting a conceptual model that identifies the key components and processes of karst aquifers. An overview for simulating the key flow and transport processes in epikarst, vadose and phreatic zones were then reviewed. Different techniques for modelling the geometrical and topological characteristics of karst conduit networks were also surveyed. When they properly account for flow and transport dynamics and accurately capture the dynamic regime transition of the studied karst system, modelling approaches can be adopted to reduce uncertainties in the management of karst groundwater. Up to now, many hydraulic behaviors and their activating thresholds that were identified by field observations and measurements have not been incorporated into physically baased karst models. Proper parameterization and efficient implementation strategies for these processes are needed to derive advanced karst models for simulating more realistic and complete hydrogeological responses. Other major issues that need to be addressed in the future are those related to development of advanced karst models integrating the complete hydraulic processes, but also the creation of realistic karst networks, the upscaling of physically baased models for catchment-scale practical applications, the study of coupled bio-hydro-chemical effects on karst aquifer evolution, linking functional and physically baased approaches. Another challenge will be the transfer of models to ungauged karst catchments on the basis of specific models linked to the typology of karst hydrological response..