Modeling the Matrix-Conduit Exchanges in Both the Epikarst and the Transmission Zone of Karst Systems

: Usual conceptual models of karst hydrodynamics highlight the important role of unsaturated subsystems in recharge repartition. However, few of them have been compared with scarce suitable physically-based numerical models. Hybrid models that couple single continuum medium with discrete features promise an improved consideration of karst speciﬁcities. Here we evaluate their capability to properly reproduce interactions between a vertical conduit and the surrounding unsaturated matrix. We simulate the response of such a conﬁguration to a single recharge event for various sets of parameters. We show the ability of hybrid models to reproduce the most signiﬁcant behaviors described in the literature, i.e., transient storage and distribution of recharge, ﬂow concentration towards conduits in the epikarst, and matrix-conduit exchanges varying in time and space. In addition to the explicit conduits, simulating variably saturated ﬂows with the Richards equation and distinguishing the epikarst and the transmission zone are key elements to reproduce most processes. The contrasts between subsystems necessary to observe desired behaviors have been quantiﬁed. They are reinforced by the varying matrix saturation that causes realistic competition between matrix and explicit conduits. The study also highlights the need to deepen knowledge of the scaled medium properties we need to know to apply such models to actual cases.


Introduction
Karst aquifers differ from other types of hydrogeological systems in their complex behavior due to strong heterogeneity. From a very tight matrix to very large conduits, including structural features at all scales, the various types of heterogeneity that constitute the karst medium induce a great diversity of interacting fluid flows and raise many questions for the hydrogeologists. Karst conduits add complexity by forming networks of preferential pathways able to drain huge, hierarchized, and possibly turbulent flows.
The study of flow processes in such complex media led numerous authors to propose conceptual models of functioning karst systems [1][2][3][4]. The usual binary classification in "conduit flow" and "diffuse flow" [5,6] was soon recognized as insufficient [4] even if both fast and slow components assume the major parts of the hydrodynamics [7]. Smart and Hobbs [8] constructed a classification including a recharge source and aquifer storage as two additional coordinates [4]. Morales et al. [9] highlighted the need to account for transient storage processes that a delayed response or a delayed infiltration can also assume [10,11]. Finally, the unsaturated zone [12][13][14] and exchanges between capacitive and conductive parts of the karst medium [4,10] appear to be some keys to these transient processes.
However, most conceptual models are mainly qualitative. They are based on field observations and lack numerical validation, likely because of scarce efficient physically-based numerical models. Few proposed conceptual models of karst unsaturated zones have been confronted with physically-based numerical models [13], rarely in 3D or with explicit conduits. Moreover, very few studies take into account unsaturated flows in karst groundwater modeling; even fewer distinguish the epikarst from the transmission zone [15][16][17][18][19][20][21]. On the other hand, when they are considered, exchanges or coupling coefficients between the fast and slow components are delicate to conceptualize [10] or to simulate [22,23] and involve crucial parameters to model calibration [23,24].
The development of numerical models capable of encompassing all the karst specificities at various scales remains an issue of interest to scientists and water managers. Different types of distributed models, more especially physically-based 3D models, have been developed [2,25,26]: equivalent porous medium approaches (EPM), discrete fracture network approaches (DFN) or discrete channel network approaches (DCN), double continuum approaches (DC), and hybrid models (HM). Double continuum models are based on two superimposed equivalent porous media representing respectively the matrix and the conduits. Usually, Darcy law applies in both media and an exchange term controls flows between the two media [19]. Although these models make it possible to distinguish the karst features from the porous matrix, they offer an unrealistic representation of the conduits and therefore a delicate parametrization of the second medium. Hybrid models offer a more realistic representation of the karst medium. They couple an equivalent porous medium with discrete features representing respectively a matrix with fractures and the networks of karst conduits. Each mesh cell is an equivalent porous medium representation of a continuum encompassing both porous rock and minor discontinuities such as small fractures and conduits that permit single medium up-scaling of petrophysical properties at the cell scale. Larger conductive discontinuities that are not compatible with up-scaling rules at the cell scale, such as main karst conduits, are considered as discrete features. Note that the cell scale should therefore largely determine the distribution of conductive objects according to whether they are supported by cells or by discrete features [27]. According to the numerical method, the edges and the faces of the mesh cells or a distinct mesh can support these discrete features. Theoretically, multiple kinds of discrete features could be considered, depending on geometry or flow physics for instance. Hybrid models hold promise for simulating actual physical processes in conduits such as turbulent flows and their interactions with the matrix continuum [23,[28][29][30]. Several authors have reported their ability to simulate the hydrodynamics of karst aquifer [31][32][33][34]. Recent tools can simplify hybrid models application to actual cases [35,36]. However, they require both data and computational capabilities that are rarely available [23,25].
This paper addresses the question of simulating flow processes in the karst unsaturated zone with hybrid models. It focuses on the interaction between explicit karst conduits and the matrix in the unsaturated zone thanks to numerical experiments. We investigated this interaction at the conduit scale in a neighboring vertical conduit. To this end, we built a hypothetical 3D hybrid model based on commonly accepted representations of the karst medium with vertical division in the epikarst, transmission zone, and saturated zone. It makes it possible to compare conceptual and simulated flow processes as a function of space and time at the conduit scale. We evaluated the ability of this hybrid model to reproduce both the percolation of water through each unsaturated subsystem and the conduit and the local exchanges between them, as compared to usual conceptual models. This raises several questions: Are simulated pressure fields and flow repartition consistent with usual representations? How do the flow repartition and the exchanges between the epikarst, transmission zone, and conduit vary as a function of space and time during a recharge event? How do varying model parameters affect flow processes and the overall simulation results? What are the model key features? What new insight do hybrid models give into the functioning of the karst unsaturated zone? What are the limitations and the ways forward of such approaches? The next section provides a literature summary of the main concepts of karst hydrodynamics and related properties. It will serve as a basis to assess the simulation results. The common properties of the modeling approach and the flow equations are also introduced.

Karst Subsystems
Hydrogeologists commonly subdivide the karst medium into four main vertical subsystems on the basis of geomorphology, hydrodynamic properties, and flow processes [1]: (i) soil and non-karst terrain, (ii) the epikarst [11,37,38], (iii) the transmission zone [39], and (iv) the saturated zone. Thus, the unsaturated zone encompasses soil, the epikarst, and the transmission zone ( Figure 1). Though this conceptual model has the advantage of being relatively generic, it remains primarily qualitative [40,41]. Some papers contributed to specify hydrodynamic behaviors for each subsystem [10,42,43]. However, most authors face difficulty in accessing the input and output of each subsystem. Particularly, the functioning of the transmission zone is generally not distinguishable from one epikarst [43] and/or one saturated zone [42].
The epikarst is the near-surface and altered karst zone [37]; its thickness ranges from 0 to approximately 20 m [37,38,44,45]. Alteration processes increase porosity and hydraulic conductivity [11,38]. Hydraulic conductivity may be higher than 10 −5 m/s and effective porosity ranges from 1% to 10% [46,47]. Due to its relatively high porosity and hydraulic conductivity compared to those of the underlying rock, the epikarst can store a large amount of water in locally perched aquifers that generally drain horizontally towards preferential flow paths [11,48,49]. Three factors affect water storage in the epikarst: (1) its thickness and lateral continuity, (2) its average porosity, and (3) the relative rates of inflow and outflow of water [45].
The transmission zone connects and transfers water from the epikarst to the saturated zone with flow essentially vertical through fractures or caves at the subsystem scale. Hydraulic conductivity for the transmission and saturated zones generally ranges between 10 −7 m/s and 10 −4 m/s, and porosity is usually less than 2% [45]. The transmission zone may also play a storage role [50,51], releasing water during the low water period [52].
Assessing the equivalent hydraulic conductivity on a flow unit of such media is not straightforward. Most of the values found in the literature (Table 1) are not related to a given measurement volume, nor the medium under consideration, which may include all or parts of the fractures and conduits. Moreover, very few studies address anisotropy even qualitatively, and even less vertical anisotropy, which largely controls the occurrence of perched and confined layers [53]. Matrix permeability is recognized to be anisotropic and higher in the dipping direction due to geological layering, with the vertical anisotropy ratio (permeability normal to the strata over permeability concordant with the strata) usually between 0.01 and 1 [54]. Alteration processes, such as fracturing, weathering, or dissolution, likely reduce or possibly reverse this vertical anisotropy of hydraulic conductivity. Due to the confinement of most of the fractures to mechanical units in sedimentary rocks, fractures are systematically wider than they are high and, except larges faults, rarely cut horizons inducing matrix vertical anisotropy [54]. Therefore, fractures increase horizontal anisotropy but affect the vertical anisotropy to a limited extent, notably in thick flow units, such as the transmission zone. It is generally assumed that karst genesis increase mainly the vertical conductivity in the unsaturated zone. However, several processes such as subaerial exposure, dissolution at the water table, or the development of karst inception features can also increase the horizontal permeability in the unsaturated zone. For instance, karst features such as lapies, exemplify the near-surface increase in vertical hydraulic conductivity. They also induce an increase in horizontal permeability and horizontal anisotropy. The net effect of these changes on the resulting vertical anisotropy is uncertain. Hence, it may be considered that: (i) the vertical karst conduits support a large part of the vertical hydraulic conductivity, (ii) these features aside, the epikarst is less anisotropic than the rest of the unsaturated zone [48]. Residual water content θ r (-) or Residual water saturation S r (-) θ r = S r = 0 [18] θ r ∈(0.01; 0.05) [15] S r = 0.05 [19] θ r = 0.171 [17] S r = 0.05 Table 1 provides a literature review of ranges for hydrodynamic parameters in karst subsystems. It serves as a basis for defining parameters of the models of this study. These parameters are also presented in Table 1.

Matrix vs. Conduit Properties
Whatever the subsystem, the karst medium is a complex heterogeneous and discontinuous spatial distribution of more or less connected entities; these entities range across a continuum from the less fractured porous matrix with high capacity and low hydraulic conductivity to karst conduits with low capacity and high hydraulic conductivity. Microporous matrix hydraulic conductivity usually ranges from 10 −12 to 10 −6 m/s (assuming that permeability of 1 Darcy is equivalent to a hydraulic conductivity of 10 −5 m/s) at a small scale [67,68], while the conduits hydraulic conductivity is greater than 10 −1 m/s [69].
Modelers usually distinguish slow or diffuse and fast or concentrated flows in all karst subsystems. The overall fast to slow flow ratio depends on karst maturity. The portion of concentrated recharge into conduits ranges from 5% for young karst to 40-58% for mature karst [13,63,65,70,71]. However, this ratio is difficult to quantify because exchanges between capacitive and conductive parts of the karst system vary as a function of space and time. Chang et al. [71] highlight the important role of the epikarst in the recharge repartition and show its impact on hybrid models simulations.

Matrix vs. Conduit Properties
Whatever the subsystem, the karst medium is a complex heterogeneous and discontinuous spatial distribution of more or less connected entities; these entities range across a continuum from the less fractured porous matrix with high capacity and low hydraulic conductivity to karst conduits with low capacity and high hydraulic conductivity. Microporous matrix hydraulic conductivity usually ranges from 10 −12 to 10 −6 m/s (assuming that permeability of 1 Darcy is equivalent to a hydraulic conductivity of 10 −5 m/s) at a small scale [67,68], while the conduits hydraulic conductivity is greater than 10 −1 m/s [69].
Modelers usually distinguish slow or diffuse and fast or concentrated flows in all karst subsystems. The overall fast to slow flow ratio depends on karst maturity. The portion of concentrated recharge into conduits ranges from 5% for young karst to 40-58% for mature karst [13,63,65,70,71]. However, this ratio is difficult to quantify because exchanges between capacitive and conductive parts of the karst system vary as a function of space and time. Chang et al. [71] highlight the important role of the epikarst in the recharge repartition and show its impact on hybrid models simulations.
It is generally considered that the conduits supply the continuum during high flow periods ( Figure  1a) whereas they drain the surrounding continuum during low flow periods ( Figure 1b) [10,11,48,72]. Such exchanges depend directly on local pressure gradients. Characteristics of conduits (e.g., diameter, hydraulic conductivity) affect also hydraulic head distribution [73][74][75][76]. However, increasing the flow capacity of the conduits above a certain threshold in saturated conditions does not affect the response of the model [74].

Modeling Approach
The following sections present the 3D hybrid model constructed to study the interactions between matrix and conduits during a recharge event. This finite element model was developed with FEFLOW 7.0 by DHI WASY (https://www.mikepoweredbydhi.com/products/feflow) [35,77].

Modeling Approach
The following sections present the 3D hybrid model constructed to study the interactions between matrix and conduits during a recharge event. This finite element model was developed with FEFLOW 7.0 by DHI WASY (https://www.mikepoweredbydhi.com/products/feflow) [35,77].

Model Description, Boundary Conditions, and Evaluation Criteria
The model represents the matrix surrounding a vertical conduit along the unsaturated zone and a part of the saturated zone (Figure 1c,d). The vertical conduit is located in the center of a 1 km 2 square so that the lateral boundaries are far enough from the conduit to have little effect on near-conduit exchanges. Hence, the presented results focus on a smaller area, up to 250 m from the conduit (e.g., Figure 2). The model is 140 m thick (from z = −40 m to z = 100 m) including the epikarst, transmission zone, and the upper part of the saturated zone. The model grid has 1405 vertex per slice and 104 slices. From the top of the model to a depth of 100 m, the layers are 1 m thick. The three deeper layers that correspond to the permanently saturated interval are thicker. They constitute a buffer zone to avoid a significant impact of the boundary conditions. The horizontal mesh is refined all around the vertical conduit.
In terms of boundary conditions, a recharge flux is applied on the top of the model (z = 100 m) whereas Dirichlet boundary conditions control the discharge at the base of the model (z = −40 m). The recharge is uniform, i.e., without a focused recharge point on discrete features, but varies as a function of time. The same single recharge event is simulated for the various parameter sets. It consists of 100 mm in two days (the third and fourth days) represented by an isosceles triangle reaching 100 mm/day in one day ( Figure 3a). This single recharge event adds to a constant recharge of 0.5 mm/day which steady-state response constitutes the initial condition of the transient simulation. This single event is realistic as observed in the Mediterranean climate for instance [78,79]. At the base of the model, constant hydraulic heads of 3 m and 2 m are applied to vertex corresponding to matrix and conduits, respectively ( Figure 1d). This difference between matrix and conduit ensures preferential flow towards the conduit at the base of the model for mean steady recharge. No-flow boundary conditions are applied to the vertical faces of the model.
The simulations results are assessed by qualitative comparison to conceptual models from the literature and quantitative inter-comparison.

Flow Equations and Model Parameters
Hybrid models make it possible to consider independent flow physics for the continuum (e.g., Darcy's law or Richards' equation) and discrete features (e.g., Hagen-Poiseuille, Darcy-Weisbach, Manning-Strickler equations, or Darcy's law). The most complex but also most realistic case would couple the Richards' equation within the continuum to the Manning-Strickler equation in the discrete features [35,70,80]. Several tests revealed that this configuration hinders convergence of the model. In this study, although failure to account for turbulent flows in discrete features may affect the overall response [29,30], Richards' equation and Darcy's law are considered for continuum and discrete features, respectively, as an acceptable compromise between numeric efficiency and realism. Flows in the matrix are assumed to obey the Richards' equation [81], which makes it possible to reproduce the variably saturated water flow: for the case when anisotropy and coordinate axes are parallel, where θ is the volumetric water content (-), t is time (s), x, y, and z are the spatial coordinates (m) (positive upwards), K(Ψ) is the unsaturated hydraulic conductivity (m/s) in the function of the pressure head Ψ (m), h is the hydraulic head (m), with h = Ψ + z, and U is the sink-source term (s −1 ). Solving the Richards' equation requires constitutive relationships for saturation as well as the relative hydraulic conductivity. Several models based on statistical pore-size distribution have been developed [82][83][84][85][86][87]. Their application and calibration mostly concern unconsolidated porous media [88]. The literature contains a few applications of the Van Genuchten model [85] to a carbonate matrix [15,18,19,21]. However, the application of this empirical model to a fractured or karstified carbonate matrix is not well documented. In fact, estimating the relationship between saturation and relative hydraulic conductivity for karst media remains a challenge [89]. Here the Van Genuchten model is therefore applied with constant and uniform parameters from the literature for all simulations (Table 1). For the unsaturated zone, water content is equal to: where θ r and θ s are residual and saturated water contents (−), respectively, and α (cm −1 ), n, and m are empirical parameters. Note that setting either residual water content θ r or residual saturation S r is equivalent because the two parameters are linked: moisture content θ equals porosity φ multiplied by saturation S. However, some confusion may exist in the literature with residual saturation values directly set as equal to residual water content values [19] ( Table 1). The relative hydraulic conductivity K r (-) in the unsaturated zone follows this relation: S e is effective saturation, generally defined as [85]: In the model, conduits are always fully conductive and follow Darcy's law. Thus, flux through the conduit is proportional to both cross-section and hydraulic conductivity of the conduit. Variations of these two parameters are therefore redundant. We thus consider the flow capacity, defined by the product of the cross-section area by hydraulic conductivity, as the key parameter for flow in conduits. The discrete features can be regarded as preferential pathways embedded in the equivalent porous medium. Discrete features and porous medium are supported by the same vertex where hydraulic heads are computed. Flow equations are solved simultaneously for both discrete features and regular mesh as if edges supporting discrete features were additional mesh cells. Exchanges between these two objects are therefore implicit [35].
We selected the values of model properties (Table 1) based on the literature. We focused on the effect of the following varying parameters: thickness (Thk), porosity (Φ), and horizontal hydraulic conductivity (K) of the epikarst (EK) and transmission zone (TZ) respectively, and flow capacity of the conduit. Values for the epikarst (Φ EK , K EK ) are distinguished from those for the transmission and saturated zone that are equal (Φ TZ−SZ , K TZ−SZ ). Assuming that the discrete feature supports most of the vertical conductivity induced by the karst and considering the contrast of anisotropy between the epikarst and transmission zone, hydraulic conductivity is set isotropic in the epikarst whereas it is set anisotropic everywhere else, with a horizontal conductivity ten times greater than the vertical, as frequently considered [54]. Note that the thickness of the epikarst (Thk EK ) does not vary during the simulation whereas the thicknesses of the transmission zone and saturated zone may vary because of flooding of the epiphreatic zone. We thus consider the initial thickness of the transmission zone as the parameter of interest (Thk TZ ). These parameters are modified around a reference simulation.
The parameter values are in the range defined in Table 1. Only one parameter is changed at a time.
Setting the flow capacity of the conduit is a tricky process. Indeed, in hybrid models, discrete elements are of interest especially for representing major conduits with high flow capacity. However, two conductive discrete features cause a matrix-restrained flow regime and act as fixed-head boundary conditions that would be more computationally efficient [74]. Thus, discrete features are of interest only if their flow capacity ensures both a high contrast with matrix characteristics and a conduit-influenced flow regime as defined by Kovács et al. [74]. Here, considering the scale of the model and some preliminary simulations, we tested values between 10 −2 and 10 1 m 3 /s −1 (Table 1 and Figure 4).

Spatio-Temporal Evolution of the Flows at the Conduit Scale
We investigated the spatial and temporal variation of the governing flow processes during the recharge event based on cross-sections of the saturation and hydraulic head. Figure 2 shows these profiles for the reference simulation at different times: at initial steady-state, at the end of the recharge event, and 77 days after the recharge event. Whatever the time, in the epikarst, water flows horizontally towards the conduit whereas the flow through the transmission zone (TZ) is mainly vertically oriented, as highlighted in the literature [1,3,48] and summarized in Figure 1a,b. The velocity field could therefore seem counter-intuitive given the anisotropy of each subsystem. However, it is mainly driven by both larger scale heterogeneity due to the contrast between subsystems and pressure at the base of the transmission zone and along the conduit.
Water 2020, 12, x FOR PEER REVIEW 8 of 20 We investigated the spatial and temporal variation of the governing flow processes during the recharge event based on cross-sections of the saturation and hydraulic head. Figure 2 shows these profiles for the reference simulation at different times: at initial steady-state, at the end of the recharge event, and 77 days after the recharge event. Whatever the time, in the epikarst, water flows horizontally towards the conduit whereas the flow through the transmission zone (TZ) is mainly vertically oriented, as highlighted in the literature [1,3,48] and summarized in Figure 1a,b. The velocity field could therefore seem counter-intuitive given the anisotropy of each subsystem. However, it is mainly driven by both larger scale heterogeneity due to the contrast between subsystems and pressure at the base of the transmission zone and along the conduit. The flow through the transmission zone is diffuse and causes the water level on the top of the saturated zone to increase up to almost 10 cm far from the conduit (Figure3d). The water level near the conduit increases also due to concentrated flow in the conduit, which slightly accentuates the piezometric dome near the conduit (Figure 3c).
Exchanges between the matrix and conduit also vary with depth and time. As expected, the conduit always drains the matrix at the top of the model. Despite a lower fixed head at the conduit boundary that could induce drainage through the whole conduit, the conduit partially supplies the matrix in the transmission zone. The depth of inversion of flux between the matrix and conduit remains relatively shallow. Even in steady-state, the conduit does not drain the matrix in the epiphreatic zone. The constant recharge probably prevents reaching the low flow as shown in Figure  1b. However, at the end of the recharge event, although the conduit is highly conductive and the boundary conditions at the bottom of the model differ between conduit and matrix, the matrix drains the conduit throughout almost the entire transmission zone. This result contrasts with the common representation of conduits draining water toward the saturated zone, probably largely influenced by results in saturated conditions [74]. The zone influenced by this flow from the conduit towards the matrix is nevertheless limited to the close vicinity of the conduit. Figure 3 confirms observations from the profiles. In the epikarst, the vertical flux through the horizontal section of the matrix decreases rapidly as a function of depth (Figure 3b) in favor of the flow transited by the conduit (Figure 3d). Consequently, the flux transited by the matrix is about forty times lower in the top of the transmission zone (Figure 3c) than in the epikarst. The transmission zone delays and flattens the matrix flux, which causes temporary storage in the medium. This is highlighted in Figure 3e which presents the flux balance of the matrix in the transmission zone, i.e., The flow through the transmission zone is diffuse and causes the water level on the top of the saturated zone to increase up to almost 10 cm far from the conduit (Figure 3d). The water level near the conduit increases also due to concentrated flow in the conduit, which slightly accentuates the piezometric dome near the conduit (Figure 3c).
Exchanges between the matrix and conduit also vary with depth and time. As expected, the conduit always drains the matrix at the top of the model. Despite a lower fixed head at the conduit boundary that could induce drainage through the whole conduit, the conduit partially supplies the matrix in the transmission zone. The depth of inversion of flux between the matrix and conduit remains relatively shallow. Even in steady-state, the conduit does not drain the matrix in the epiphreatic zone. The constant recharge probably prevents reaching the low flow as shown in Figure 1b. However, at the end of the recharge event, although the conduit is highly conductive and the boundary conditions at the bottom of the model differ between conduit and matrix, the matrix drains the conduit throughout almost the entire transmission zone. This result contrasts with the common representation of conduits draining water toward the saturated zone, probably largely influenced by results in saturated conditions [74]. The zone influenced by this flow from the conduit towards the matrix is nevertheless limited to the close vicinity of the conduit.  (Figure 3d). Consequently, the flux transited by the matrix is about forty times lower in the top of the transmission zone (Figure 3c) than in the epikarst. The transmission zone delays and flattens the matrix flux, which causes temporary storage in the medium. This is highlighted in Figure 3e which presents the flux balance of the matrix in the transmission zone, i.e., the difference between the flux into and out of the matrix in the transmission zone, as a function of time. Exchanges between the matrix and conduits are limited or almost balanced in the transmission zone. Indeed, whatever the time, the flux transited by the conduit at the bottom (5 m elevation) of the transmission zone is slightly smaller than the one at the top (80 m elevation) of this zone (Figure 3c) and the flux balance in the conduit (the difference between the flux into and out of the conduit in the transmission zone) is always smaller than 0.21 L/s −1 (Figure 3f). Here, water from the conduit wet the unsaturated medium surrounding the conduit. However, the total flux from the conduit to the transmission zone remains below 3% of the flux entering the transmission zone by the conduit. As observed on the profiles, the influence on flows in the transmission zone remains limited to the area close to the conduit.
At the base of the epikarst, far from the conduit, saturation increases rapidly during the recharge event and decreases slowly after it ( Figure 2). The vertical conduit finally drains the water temporarily stored there. This configuration does not produce a locally perched aquifer, as referred to in the literature. Nevertheless, these results highlight several factors that may favor the occurrence of locally saturated zones in the epikarst: the distance from the conduit and contrasts of hydraulic conductivity. More important (repeated, longer, or more intense) recharge events and topographic variations of the interface between the epikarst and the transmission zone likely favor these phenomena too. The impact of varying parameters on saturation at the bottom of the epikarst is discussed hereafter.

Impact of Parameter Variation at the Conduit Scale
The effect of the extreme values of the flow capacity of the conduit is limited. Obviously, a too-small value of flow capacity (here, 10 −2 m 3 /s −1 ) limits the preferential flow through the conduit (Figure 4). Conversely, increasing the flow capacity of the conduits above a threshold value (here, 1 m 3 /s −1 ) does not have a greater impact on the response of the model, as referred to by Kovács et al. [74] for saturated conditions. We have therefore selected an intermediate value (10 −1 m3/s −1 ) for the reference model, to simulate a highly conductive conduit that ensures preferential flow but whose effect differs from that of a fixed-head boundary condition. Figures 5-8 show profiles of the saturation and hydraulic head for different varying parameters. Each profile is taken at the time step for which the flux at the bottom output of the model is maximum. Figure 5 shows four profiles of the saturation and hydraulic head for different values of K EK (from 10 −5 to 10 −1 m/s) while K TZ−SZ (10 −5 m/s) and other parameters remain constant. Conversely, Figure 6 shows three profiles of the saturation and hydraulic head for different values of K TZ−SZ (from 10 −7 to 10 −3 m/s) while K EK and other parameters remain constant (10 −2 m/s). The contrast of hydraulic conductivity between the epikarst and transmission zone mainly controls the flow direction in the epikarst during the recharge event. Other things being equal, a ratio of hydraulic conductivity higher than ten (one hundred for the ratio of vertical conductivity due to the difference of anisotropy) ensures horizontal preferential flow through the epikarst toward the conduit while saturation monotonically increases as a function of depth (Figure 5a versus Figure 8b,c). A ratio greater than about 10 4 is necessary to distinguish clearly specific behaviors in the epikarst as compared to the transmission zone: local saturation is higher at the bottom of the epikarst, sharp variation of flow direction between both zones, and the conduit is supplying the matrix in the transmission zone (Figures 5d and 6a).  The effect of the extreme values of the flow capacity of the conduit is limited. Obviously, a toosmall value of flow capacity (here, 10 −2 m 3 /s −1 ) limits the preferential flow through the conduit ( Figure  4). Conversely, increasing the flow capacity of the conduits above a threshold value (here, 1 m 3 /s −1 ) does not have a greater impact on the response of the model, as referred to by Kovács et al. [74] for saturated conditions. We have therefore selected an intermediate value (10 −1 m3/s −1 ) for the reference model, to simulate a highly conductive conduit that ensures preferential flow but whose effect differs from that of a fixed-head boundary condition.   Figure 5 shows four profiles of the saturation and hydraulic head for different values of KEK (from 10 −5 to 10 −1 m/s) while KTZ−SZ (10 −5 m/s) and other parameters remain constant. Conversely, Figure 6 shows three profiles of the saturation and hydraulic head for different values of KTZ−SZ (from 10 −7 to 10 −3 m/s) while KEK and other parameters remain constant (10 −2 m/s). The contrast of hydraulic conductivity between the epikarst and transmission zone mainly controls the flow direction in the epikarst during the recharge event. Other things being equal, a ratio of hydraulic conductivity higher than ten (one hundred for the ratio of vertical conductivity due to the difference of anisotropy) ensures horizontal preferential flow through the epikarst toward the conduit while saturation monotonically increases as a function of depth (Figure 5a versus Figures 8b,c and 9b,c). A ratio greater than about 10 4 is necessary to distinguish clearly specific behaviors in the epikarst as compared to the Varying porosity in the epikarst (Figure 7) also affects the flow and saturation repartition in the entire medium. Although it is difficult to compare the logarithmic variation of hydraulic conductivity with the linear variation of porosity, the impact of the latter is greater in the range of realistic values for both parameters. It is presumably due to the combined effects of porosity on both saturation and effective conductivity, in accordance with the application of the Richards' equation and the Van Genuchten empirical model. Indeed, for the same amount of incoming water, the lower the porosity values, the higher the saturation and the effective conductivity. Lower porosity values thus yield lower local effective storage and transit time.
Various values of thickness were also tested for the epikarst and the transmission zone. Figure 8 shows profiles of the saturation and hydraulic head for three different values of thickness of the epikarst. The thickness of the entire system remains constant and the thickness of the transmission zone varies consequently. Note the specific case of a null thickness corresponding to the absence of the epikarst (Figure 8a). The three results are very different. For a sufficient thickness, here above 10 m, the profiles have all the features mentioned above, notably horizontal flow towards the conduit in the epikarst, the supply of the matrix through the conduit at depth, and a higher saturation far from the conduit at the base of the epikarst. A decreasing thickness of the epikarst causes decreasing storage capacity and, for the same recharge event, the saturation and hydraulic head reach higher values with higher gradients. The flow features (horizontal flow towards the conduit in the epikarst, the supply of the matrix through the conduit at depth) are preserved even for a thin epikarst (Figure 8b) because the interface between the epikarst and the transmission zone acts as a hydraulic barrier. Indeed, the saturation contrast created by the recharge event enhances the conductivity contrast between the two subsystems. Thus, the top of the transmission zone remains poorly saturated, therefore not very conductive, while preferential flows drain water towards the conduit in the more conductive epikarst. However, the shape of the saturation profiles varies with the epikarst thickness. The formation of a zone with higher saturation far from the conduit at the base of the epikarst requires a sufficient epikarst thickness. Without the epikarst (Figure 8a), the resulting profiles are very different, with higher saturation values at the top of the transmission zone, a monotonous variation of saturation, and the repartition of flows towards the conduits all along the transmission zone. Thus, notably due to unsaturated conditions and whatever its thickness, the epikarst appears to be a key feature for concentrating water towards conduits in addition to geomorphology and topography which also promote such processes in real karst systems.
Water 2020, 12, x FOR PEER REVIEW 12 of 20 transmission zone: local saturation is higher at the bottom of the epikarst, sharp variation of flow direction between both zones, and the conduit is supplying the matrix in the transmission zone (Figures 5d and 6a).   transmission zone: local saturation is higher at the bottom of the epikarst, sharp variation of flow direction between both zones, and the conduit is supplying the matrix in the transmission zone (Figures 5d and 6a).   with the linear variation of porosity, the impact of the latter is greater in the range of realistic values for both parameters. It is presumably due to the combined effects of porosity on both saturation and effective conductivity, in accordance with the application of the Richards' equation and the Van Genuchten empirical model. Indeed, for the same amount of incoming water, the lower the porosity values, the higher the saturation and the effective conductivity. Lower porosity values thus yield lower local effective storage and transit time. Various values of thickness were also tested for the epikarst and the transmission zone. Figure  8 shows profiles of the saturation and hydraulic head for three different values of thickness of the epikarst. The thickness of the entire system remains constant and the thickness of the transmission zone varies consequently. Note the specific case of a null thickness corresponding to the absence of the epikarst (Figure 8a). The three results are very different. For a sufficient thickness, here above 10 m, the profiles have all the features mentioned above, notably horizontal flow towards the conduit in the epikarst, the supply of the matrix through the conduit at depth, and a higher saturation far from the conduit at the base of the epikarst. A decreasing thickness of the epikarst causes decreasing storage capacity and, for the same recharge event, the saturation and hydraulic head reach higher values with higher gradients. The flow features (horizontal flow towards the conduit in the epikarst, the supply of the matrix through the conduit at depth) are preserved even for a thin epikarst ( Figure  8b) because the interface between the epikarst and the transmission zone acts as a hydraulic barrier. Indeed, the saturation contrast created by the recharge event enhances the conductivity contrast between the two subsystems. Thus, the top of the transmission zone remains poorly saturated, therefore not very conductive, while preferential flows drain water towards the conduit in the more conductive epikarst. However, the shape of the saturation profiles varies with the epikarst thickness. The formation of a zone with higher saturation far from the conduit at the base of the epikarst requires a sufficient epikarst thickness. Without the epikarst (Figure 8a), the resulting profiles are very different, with higher saturation values at the top of the transmission zone, a monotonous variation of saturation, and the repartition of flows towards the conduits all along the transmission zone. Thus, notably due to unsaturated conditions and whatever its thickness, the epikarst appears to be a key feature for concentrating water towards conduits in addition to geomorphology and topography which also promote such processes in real karst systems. Varying the thickness of the transmission zone provides similar profiles of the hydraulic head (not presented here), with lower heads for lesser thicknesses. Higher thickness of the transmission zone causes lower saturation and consequently lower effective conductivity at the top of this zone, in accordance with Richards' equation. Therefore, the hydraulic conductivity contrast at the interface between the epikarst and transmission zone is higher, which favors the development of a local maximum saturation at the base of the epikarst. The drainage radius of the conduit and flow anisotropy also increases as the thickness of the transmission zone increases.
Finally, the transmission zone largely affects the distribution of pathways and transit times between the matrix and conduits. On the one hand, by contrast with the epikarst properties, it could act or not as a barrier that reinforces horizontal drainage towards conduits in the epikarst. On the other hand, the transmission zone directly competes with the vertical conduit to flow water from the epikarst to the saturated zone.

Discussion
We performed numerical experiments to evaluate the ability and interest of hybrid models to represent flow processes in subsystems of the karst unsaturated zone and the associated matrixconduits exchanges. We built the model with the aim of reproducing commonly accepted conceptual models while setting parameter values consistent with the literature. The assessment of model results is based on qualitative comparisons of cross-sections of simulated saturation and hydraulic head with conceptual representations from the literature illustrated in Figure 1. Our results show that hybrid Varying the thickness of the transmission zone provides similar profiles of the hydraulic head (not presented here), with lower heads for lesser thicknesses. Higher thickness of the transmission zone causes lower saturation and consequently lower effective conductivity at the top of this zone, in accordance with Richards' equation. Therefore, the hydraulic conductivity contrast at the interface between the epikarst and transmission zone is higher, which favors the development of a local maximum saturation at the base of the epikarst. The drainage radius of the conduit and flow anisotropy also increases as the thickness of the transmission zone increases.
Finally, the transmission zone largely affects the distribution of pathways and transit times between the matrix and conduits. On the one hand, by contrast with the epikarst properties, it could act or not as a barrier that reinforces horizontal drainage towards conduits in the epikarst. On the other hand, the transmission zone directly competes with the vertical conduit to flow water from the epikarst to the saturated zone.

Discussion
We performed numerical experiments to evaluate the ability and interest of hybrid models to represent flow processes in subsystems of the karst unsaturated zone and the associated matrix-conduits exchanges. We built the model with the aim of reproducing commonly accepted conceptual models while setting parameter values consistent with the literature. The assessment of model results is based on qualitative comparisons of cross-sections of simulated saturation and hydraulic head with conceptual representations from the literature illustrated in Figure 1. Our results show that hybrid models are capable of reproducing most accepted features of karst hydrodynamics in an unsaturated zone: storage and distribution of recharge in the epikarst, flow concentration towards conduits in the epikarst, intensity and direction of exchanges between the matrix and conduit that vary in time and space. The consideration of an epikarst layer and the contrast of effective conductivity between the epikarst and transmission zone are key elements for enhanced flow simulation in the karst unsaturated zone. Moreover, this model gives some insights to understand and quantify the effect on flow processes of parameter variation in each subsystem. This opens the door to testing reservoir-scale applications including complex networks of karst conduits [90].
Using Richards' equation and the Van Genuchten model to consider variably saturated flow also appears to be very useful in reproducing these processes, including saturation that varies in time and space. Several tests not presented here confirm that saturation front displacement and perched aquifers are likely impossible to simulate with simplified representations of unsaturated flows that use pseudo-saturation [91], such as those available in FEFLOW [35]. Applying the Richards' equation to the continuum has nevertheless led us to simplify the representation of flows in the conduits without simulating turbulent flows. If the consideration of turbulent flow is needed for correctly modeling large-scale saturated aquifers with small or rough conduits [30], the configuration and the subsequent needs are completely different here: a variably saturated medium with a major conduit is considered at mesoscale. The consistency between the simulations results and both observations and concepts issued from real karst systems, as well as the wide spectrum of responses obtained, lead us to consider that the major features have probably been identified and that hybrid models allow realistic simulations of unsaturated flow processes near karst conduits. However, addressing the issue of considering turbulent flow in the conduit together with the variably saturated flow in the continuum and, in general, enhancing and upscaling physically-based modeling of flows in both the matrix and conduits of the unsaturated zone, still poses numerical challenges.
On the other hand, parameters of the Van Genuchten model have not been varied due to a lack of information in the literature. The setting of the constitutive relationships for saturation as well as the relative hydraulic conductivity, as the Van Genuchten model, applied to fractured and karstified carbonates is a broad issue whatever the scale of interest. Further progress in this area would require a dedicated study. More generally, setting the parameters of hybrid models appears to be a tricky task. Beyond the usual problems of investigation techniques, data acquisition, and up-scaling largely addressed in the literature [2,25], improving the physics represented in the models has a corollary to define and assess their parameters correctly. The continuum and the discrete elements of hybrid models each support part of the characteristics of the karst medium that must, therefore, be distributed [27]. Some questions remain and are probably case-dependent: what exactly do we represent with discrete elements, geometric objects, or physical processes and how to up-scale them? In the case of the hydraulic conductivity field, for example, it is not only difficult to attribute to each entity a conductivity value that honors the various constraints including suitability for the represented object or process but also to evaluate related parameters such as anisotropy. Additionally, in this hypothetic case, the problem of intra-subsystem heterogeneity has not been addressed. Thus, the availability of models capable of finely representing the physical processes at various scales transforms the paradigm of succeeding in modeling field observations in the question of characterizing the properties of the medium studied for modeling purposes.

Conclusions
In this paper, we examined some usual concepts of karst hydrodynamics using a hybrid model. We focused this work on interactions between a fully conductive conduit and the surrounding unsaturated matrix, which have rarely been addressed. We explored a range of parameters from the literature to evaluate the ability of such models to simulate observed processes, to quantify how varying parameters affect these processes, and to identify the key features and parameters of the model. Especially, we focused on the distinction between the epikarst and the transmission zone.
The hybrid model was able to reproduce conceptual representations of karst hydrodynamics at the conduit scale with ranges of parameters commonly accepted in the literature: transient storage and distribution of recharge, flow concentration towards conduits in the epikarst, and matrix-conduit exchanges, with intensity and direction that vary in time and space. Considering variably saturated flow appears to be useful for reproducing most of the processes referenced in the literature. Indeed, saturation varying as a function of time in the matrix also caused a realistic variability in the competition between the matrix and explicit conduits. The study also raised the need to deepen knowledge on the saturation-conductivity relationships in karst media and to develop hybrid models effectively coupling the modeling of unsaturated flows in the matrix to the modeling of turbulent flows in the discrete elements. The recent work opens interesting avenues [23]. Finally, the remaining question is more about the estimation of the scaled properties we need to know to simulate physical processes correctly than how to simulate them.
The various numerical simulations confirmed that the usual distinction between the epikarst and the transmission zone is important for controlling the near-surface flow distribution. The exchanges between the epikarst and the conduit, presented in Section 2.1 and summarized in Figure 1a,b, were properly reproduced and highlighted the important role of the epikarst in distributing the recharge. Our study confirmed that the presence of the epikarst affects not only the distribution of recharge in this zone but also the conduit-matrix interactions in the entire unsaturated zone, as observed by Trček [92]. Moreover, the variation of the effective conductivity related to the variation of saturation in the medium reinforced the contrast between the epikarst and the transmission zone during a recharge event and affected the exchange between the matrix and the conduit in close proximity to the conduit.