The Effects of Exchange Flow on the Karst Spring Hydrograph under the Different Flow Regimes: A Synthetic Modeling Approach

In this study, a synthetic modeling approach is proposed to quantify the effect of the amount and direction of the exchange flow on the karstic spring discharge fluctuations under different hydrologic conditions corresponding to high and low flow conditions. We hypothesis that the spring discharge fluctuations constitute a valuable proxy to understand the internal processes of the karst system. An ensemble of spring hydrographs was synthetically produced to highlight the effect of exchange flow by exploring the plausible range of variability of coefficients of exchange flow, conduit diameter, and matrix hydraulic conductivity. Moreover, the change of the rate of point recharge through the karst conduit allows for the quantifying of the sensibility of the spring hydrograph to the directions of exchange flow. We show that increasing the point recharge lies to a remarkable linear recession coefficient (β) as an indication of the conduit flow regime. However, a reduction in and/or lack of the point recharge caused the recession coefficient to change to exponential (α) due to the dominant effect of the matrix restrained flow regime and/or conduit-influenced flow regime. The simulations highlight that the exchange flow process from the conduit to the matrix occurred in a short period and over a restricted part of the conduit flow regime (CFR). Conversely, the exchange flow dumped from the matrix to the conduit occurs as a long-term process. A conceptual model is introduced to compare spring hydrographs’ characteristics (i.e., the peak discharge, the volume of baseflow, and the slope of the recession curve) under the various flow conditions with the directions of the exchange flow between the conduit and the matrix.


Introduction
Karst aquifers can be considered as a "dual-porosity" or "dual-permeability" media ( Figure 1) that include (1) matrix porosity comprising intergranular pores and the small joints with high storage capacity and low velocity, and (2) conduit porosity comprising enlarged, dissolved, and unclogged fractures with low storage capacity and high flow velocity [1]. Karst aquifers are characterized by a strong physical heterogeneity and may present strong anisotropy in terms of hydrodynamic properties such as hydraulic conductivity and specific storage due to a complex and three-dimensional dissolution development [2].
These conditions imply an inherent duality in all processes governing the hydrological response of the karst system, including (1) the infiltration and recharge processes which range from diffusive and slow recharge into the matrix to rapid and concentrated recharge into the conduit through features such as sinkholes, sinking streams, and ponor caves,  [16], Hartmann et al. [3], Binet et al. [33]).

Theoretical Background
To quantify the exchange flow temporal fluctuations, the linear Darcy flow equation was conventionally applied. Barenblatt et al. [7] considered two assumptions for a linear, slow, and steady exchange flow between pores and fissures: (1) the inertialess flow under Darcy's law through fissures with small size and velocity, and (2) the stationary flow with the smooth changes of hydraulic gradient between pores and fissures.
Supplementary Table S1 summarizes a list of different approaches that have been used to quantify the exchange flow. The conventional equation for calculating the exchange flow in the karst aquifers is described by Bauer et al. [15]:  [16], Hartmann et al. [3], Binet et al. [33]).

Theoretical Background
To quantify the exchange flow temporal fluctuations, the linear Darcy flow equation was conventionally applied. Barenblatt et al. [7] considered two assumptions for a linear, slow, and steady exchange flow between pores and fissures: (1) the inertialess flow under Darcy's law through fissures with small size and velocity, and (2) the stationary flow with the smooth changes of hydraulic gradient between pores and fissures.
Supplementary Table S1 summarizes a list of different approaches that have been used to quantify the exchange flow. The conventional equation for calculating the exchange flow in the karst aquifers is described by Bauer et al. [15]: where Q ex is the exchange flow rate (L 3 T −1 ), α ex is the lumped exchange coefficient (L 2 T −1 ), h c and h F are the hydraulic head in the conduit and fissured matrix (L), respectively. Although Equation (1) was suggested for the steady-state condition, the hydraulic head difference between the matrix and the conduit can be considered as time-dependent ( [16,24,34]) and spatially distributed [17]. The exchange flow tends to be negligible if the head difference between the matrix and conduits tends to zero [33]. The amount of the exchange flow depends on the exchange coefficient and the difference of hydraulic head between the matrix and the conduit (Equation (1)). Estimating the exchange coefficient between two compartments with different hydraulic properties constitutes a challenge, even if the hydraulic head can be measured in both compartments. The exchange coefficient constitutes the conceptual representation of the lumped hydraulic and geometrical properties of the karst system including the conduit and matrix [15,21,23]. Equation (2) constitutes a classical approach to calculate the exchange coefficient [15]: where ∆L is the distance traveled by flow (L), σ is a parameter that depends on conduit geometry and defines as inverse fissure spacing (L −1 ), A F is the exchange surface between the conduit and the fissured matrix (L 2 ), and K F is the hydraulic conductivity of the fissured matrix (L T −1 ). The first parameter on the right-hand side of Equation (2), σ, can be considered as a criterion of the specific surface of the fissures, (i.e., the surface of the fissures per unit of volume of the rock [7]. The second parameter, A F , represents the surface of the interaction between the matrix and conduit. A F was more specifically defined by Shoemaker et al. [35] and Peterson and Wicks [14]. In another study, Reimann et al. [20] used an exchange perimeter (L) instead of the exchange surface. The third parameter, K F , represents the matrix hydraulic conductivity which linearly influences the exchange flow [5]. A wide range of the matrix hydraulic conductivity of 10 −4 to 10 −7 m/s has been used to model karst aquifers hydrological response [19,20,33,36,37].
Moreover, Bauer et al. [15] estimated a lumped parameter, α ex , instead of estimating each of the parameters on the right-hand side of Equation (2) (i.e., σA F K F ), reducing the uncertainty of Equation (2) because of the unknown geometry of the block or the fracture. They concluded a constant exchange coefficient of 1 × 10 −7 m 2 /s for the simulation of the synthetic karst aquifer by sensitivity analysis. Malenica et al. [38] calibrated values of 0.18 and 0.23 m 2 /s for the exchange coefficient based on the different experimental setups and the change of properties of conduits such as position, saturation percentage, roughness coefficient, and hydraulic head in the upstream and downstream. The influence of the conduit geometry and network on the exchange flow was examined by Saller et al. [21] via sensitivity analysis methods. A wide range from 1 × 10 −8 [36] to 25 m 2 /s [13,37] has been proposed for the lumped exchange coefficient. However, Dufoyer et al. [19] assigned a narrow range of 0.1 to 1.5 (unitless) to the exchange coefficient for simulating a karst aquifer based on the numerical double porosity model called MARTHE. Moreover, De Rooij [39] suggested a well-index in the DisCo model instead of the lumped exchange coefficient that is defined in the MODFLOW-2005 CFPM1 code. The dependence of the well-index on the spatial discretization constitutes the main difference with the lumped exchange coefficient.

. Temporal and Spatial Variations of the Exchange Flow
According to Equation (1), the direction and amount of the exchange flow are controlled by the hydraulic gradient between the conduit and the matrix [40]. Since the hydraulic gradient between the matrix and conduit is variable over time due to recharge events and/or at different zones of karst aquifers (i.e., spatial heterogeneity), many researchers have considered the exchange flow as time and space dependent.
Zhang et al. [17] used a conceptual model to simulate bi-directional exchange flow and solute transport between the matrix and the conduit network in the karst catchment with an area of 73.5 km 2 in southwest China. This conceptual model combines the dual reservoir runoff model with the ARMA model (Auto-Regressive and Moving Average Model). They concluded that during the dry season, as the water level decreases in conduits more rapidly than that in the matrix, the storage water in the matrix mostly drains into conduits as the baseflow which is indicated by high values of EC. In contrast, in the wet season, especially during periods of recharge, conduits are rapidly recharged by the infiltrated water as well as the drop of the EC values is observed in the karst spring. As a result, the water level in the conduits is higher than the water level in the matrix. In this case, the water is temporarily stored in the conduits and then gradually transferred to the matrix. Zhang et al. [17] suggested that the spatial variation of characteristics of karst aquifers, such as hydraulic gradient between the matrix and conduits, the storage capacity, and the contact area from the upstream to downstream, lead to spatial variations in the exchange flow.
In a comprehensive study, the spatial distribution of the exchange flow from the recharge zone to the discharge zone in a karst aquifer with the allogeneic recharge was investigated by Binet et al. [33] in Val d'Orléans, France ( Figure 1).
They found three zones for the spatial distribution of exchange flow in the studied karst aquifer: (1) recharge zone in the upstream, where the direction of the exchange flow is from the matrix to the conduit, (2) the intermediate zone, where there is no exchange flow, and (3) the discharge zone in the downstream, where the exchange flow comes from the conduit to the matrix. A similar distinction was performed to study the spatial distribution of exchange flow in the unsaturated zone by Soglio et al. [41]. Moreover, the spatial and temporal variations of the exchange flow can be affected by the location and type of wells (pumping or injection) [42].

The Effect of the Exchange Flow and Its Direction on the Spring Hydrograph
A spring hydrograph illustrates the temporal variations of the spring discharge reflecting the hydraulic processes in a karst aquifer [43]. Typically, a spring hydrograph is characterized by the peak discharge, the rising and the falling limbs, and two inflection points ( Figure 2). The inflection points of the rising limb and the falling limb represent the maximum infiltration state [44] and the end of an infiltration event [43], respectively. The inflection point of the falling limb (i.e., the recession limb) is determined as the intersection of the early steep slope segment (i.e., flood recession) and the late smooth slope segment (i.e., baseflow recession).
Under various hydrological conditions, three flow regimes may influence the falling limb including (1) matrix restrained flow regime (MRFR), when conduit network is as fix-head boundary conditions, and the drainage process is controlled by the matrix alone, (2) conduit-influenced flow regime (CIFR), when the hydraulic gradient in the conduits is not negligible during the recession process and the recession process depends on the hydraulic parameters matrix and conduit network [37,45], and (3) conduit flow regime (CFR), when the hydraulic head in the conduit is greater than the hydraulic head in the matrix [16]. Generally, the flow regime in the karst aquifers is influenced by the point and diffuse recharge events and the contribution of the matrix and/or conduits through-flow system. In particular, the recession limb of the spring hydrograph constitutes a linear and non-linear form in response to large and small recharge events, respectively [46]. Bailly-Comte et al. [16] showed that a linear decrease in the karst spring's discharge can Water 2021, 13, 1189 7 of 25 be assigned to the conduit flow regime (CFR) characterized by the low permeability of the matrix and concentrated recharge to the karst conduit's network. The slope of the recession curve in the conduit flow regime is denoted by β and can be calculated based on Equation (4) [16]: where Q CFR (t) is the discharge at time t (L 3 T −1 ), Q max is the maximum discharge (L 3 T −1 ), and β (L 3 T −2 ) is the linear decrease in the recession curve slope.
Under various hydrological conditions, three flow regimes may influence the falling limb including (1) matrix restrained flow regime (MRFR), when conduit network is as fixhead boundary conditions, and the drainage process is controlled by the matrix alone, (2) conduit-influenced flow regime (CIFR), when the hydraulic gradient in the conduits is not negligible during the recession process and the recession process depends on the hydraulic parameters matrix and conduit network [37,45], and (3) conduit flow regime (CFR), when the hydraulic head in the conduit is greater than the hydraulic head in the matrix [16]. Generally, the flow regime in the karst aquifers is influenced by the point and diffuse recharge events and the contribution of the matrix and/or conduits through-flow system. In particular, the recession limb of the spring hydrograph constitutes a linear and nonlinear form in response to large and small recharge events, respectively [46]. Bailly-Comte et al. [16] showed that a linear decrease in the karst spring's discharge can be assigned to the conduit flow regime (CFR) characterized by the low permeability of the matrix and concentrated recharge to the karst conduit's network. The slope of the recession curve in the conduit flow regime is denoted by β and can be calculated based on Equation (4) [16]: where QCFR (t) is the discharge at time t (L T ), Qmax is the maximum discharge (L T ), and β (L T ) is the linear decrease in the recession curve slope.
Considering concentrated recharge events, the inflection point of the falling limb has been introduced as an indication of the completion of the infiltration and recharge processes and the CFR [16,43]. However, considering diffuse recharge events, the conduit network constitutes fixed-head boundary conditions and flow conditions correspond to a matrix flow regime, with the inflection point of the falling limb only indicating the fulfillment of the recharge and infiltration processes. After completion of the recharge (i.e., beyond the inflection point of the falling limb), the baseflow exponentially decreases over time and can be calculated based on Equation (5) [47]: Considering concentrated recharge events, the inflection point of the falling limb has been introduced as an indication of the completion of the infiltration and recharge processes and the CFR [16,43]. However, considering diffuse recharge events, the conduit network constitutes fixed-head boundary conditions and flow conditions correspond to a matrix flow regime, with the inflection point of the falling limb only indicating the fulfillment of the recharge and infiltration processes. After completion of the recharge (i.e., beyond the inflection point of the falling limb), the baseflow exponentially decreases over time and can be calculated based on Equation (5) [47]: where Q 0 and Q (t) are the discharges at the initial time and time t (L 3 T −1 ), respectively, and α (T −1 ) is the recession coefficient.
Since the baseflow of the recession curve (beyond the inflection point of the falling limb) is less influenced by the recharge events, it can be considered as representative of the general characteristics of karst aquifers. Therefore, the baseflow recession can be influenced by MRFR or CIFR [43,45].
Although the recession curve of a homogeneous matrix can be decomposed into several exponential components (α), these components should not be considered as discharge from different parts of the aquifer with different hydraulic conductivity due to the complexity of diffuse flow through the matrix [43]. Fiorillo [48] showed that the variation of the recession coefficient can be explained by the variation of the discharge area without any change in permeability of the material. High values of α can be associated by a small catchment area or by drainage of water-filled shafts.
Moreover, the direction of the exchange flow can temporally vary due to the establishment of the MRFR/CIFR or CFR conditions [16]. Although, the general direction of the exchange flow is from the matrix to conduits (i.e., the exponential slope of the recession curve), a high-pressure head in conduits is necessary to establish an exchange flow from conduits to the matrix (i.e., the linear slope of the recession curve). This high-pressure head in conduits has been simulated by increasing the concentrated recharge through conduits as suggested by Malenica et al. [38].
In addition to the direction of exchange flow, the different values of the exchange coefficient, as well as the amount of exchange flow, influence the characteristics of the hydrograph. The larger exchange coefficient allows for the easy transfer of a significant amount of water from conduits to the matrix in response to a significant rainfall event. Also, the slope of the recession curve reduces and the baseflow becomes more important when reversion of the exchange flow direction occurs from the matrix to the conduits during a baseflow period [5,18]. Therefore, the small exchange coefficients in comparison to the large exchange coefficients produce a larger peak discharge and a steeper slope of the recession segment under point recharge conditions [18,19]. The effect of the exchange flow on the discharge volume of the spring depends on the matrix hydraulic conductivity. The large matrix hydraulic conductivity allows the point recharged water to flow easily from the conduit toward the matrix. Accordingly, the quick-flow tends to decrease. In contrast, in karst aquifers with low matrix hydraulic conductivity, the quick-flow volume is affected by the characteristics of the conduits [16]. Peterson and Wicks [14] estimated that the volume of exchange flow in the Devil's Icebox-Connor's Cave System (central Missouri, USA) changed from 0.34 to 6.85 × 10 −5 m 3 if the matrix hydraulic conductivity decreased from 8.54 × 10 −6 to 1 × 10 −11 m/s, respectively.
The transition from pressurized to free-surface flow or vice versa in the conduits caused the changes in spring discharge. Although the exchange flow rate from the conduit to the matrix increased as a result of a quick increase in the water pressure in the filled conduits than partly filled conduits during rainfall, the spring discharge is less affected by the exchange flow in conduits with the free-surface flow [20].

Methodology
We hypothesize that the matrix-conduit exchange flow rate depends on the lumped exchange coefficient and the hydraulic head difference between the matrix and the conduit (Equation (1)). This study attempted to evaluate the effects of these parameters on the rate and direction of the exchange flow in a karst system using a synthetic modeling approach with MODFLOW-2005 Conduit Flow Process (CFP) developed by Shoemaker et al. [35]. The CFP package [35] allows groundwater flow simulation in the matrix based on Darcy's law and flow simulation in the conduit based on the linear and nonlinear flow regime developed by the U.S. Geological Survey, (USGS) [30].
The proposed methodology consists of (1) building a synthetic karst aquifer to simulate spring hydrographs for a physical realist range of values for the exchange coefficient, the matrix hydraulic conductivity, and the conduit diameter, and (2) assessing the effect of variations in the parameter of interest on the exchange flow in terms of rate and flow direction.
The characteristics of the spring hydrographs, including peak discharge, peak time, the slope of the recession limb, and volume of baseflow (the volume of spring discharge from inflection point to the end of simulation time) were calculated and interpreted based on the rate and direction of the exchange flow.
Three-dimensional laminar flow in the matrix is solved based on Equation ( where K xx , K yy , K zz are the hydraulic conductivity values in the x, y, z directions, respectively (L T −1 ), h is the hydraulic head [L], W is the volumetric flux per unit volume (T −1 ), S s is the specific storage (L −1 ), and t is time (T). The conduit network within the matrix porous media is defined as pipes characterized by elevation, diameter, roughness, and exchange coefficient. The two ends of the pipes are connected to the nodes at the center of each matrix cell [35]. Depending on the Reynolds number, the flow in the pipes can be either turbulent or laminar flow.
Both laminar and turbulent pipe flow can be approximated with the Hagen-Poiseuille (Equation (7)) and the Darcy-Weisbach equations (Equation (8)), respectively [35]: where d is the pipe diameter (L), ∆h ∆l is the hydraulic head gradient along the pipe (unitless), ρ is the density of the water (ML −3 ), µ is the viscosity of the water (ML −1 T −1 ), g is the gravitational acceleration constant (L T −2 ), and τ is tortuosity (unitless).
where Q is the volumetric flow rate (L 3 T −1 ), ∆h is the hydraulic head losses along the pipe (L), A is the cross-sectional area perpendicular to flow (L 2 ), f is the friction factor (unitless) and is quantified based on the empirical Colebrook-White formula [49,50] in CFPM1.
where k c (L) represents the mean roughness height of the conduit wall micro-topography. R e is Reynolds number (unitless) and d is the pipe diameter (L). Using the CFPM1 package, the equation to estimate discharge in the conduit is solved independently of the matrix equation [35]. However, the relationship between the matrix and the pipe network is shown by the influence of exchange flow as an inflow or outflow component in the water balance as well as the distribution of the hydraulic head in each of the matrix and pipe networks. For any pipe nodes, conservation volumetric flow is calculated according to Kirchhoff's law (Equation (10)).
where Q ip , Q ex , Q R , and Q s represent pipe flows, exchange flow, direct recharge, and changes in storage, respectively [35]. Briefly, the volumetric exchanges (Q ex ) are defined as the total volumetric exchange Q tex (L 3 T −1 ) as follows: where in is an integer subscript for each node, mxnode is the maximum number of nodes in the model, ip is an integer subscript for each pipe that connects to node in, np is the number of pipes connected to node in. Individual volumetric exchanges (Q ex ) are computed using the Darcy formulation assuming that exchange flow is not turbulent.
where α ex(j,i,k) is the exchange coefficient (L 2 T −1 ), h in is the head (L) at node in computed by the CFP, h j,i,k is the head (L) in the encompassing matrix cell [35]. The exchange coefficient can be determined by two options in CFP. The first option allows the user to assign the exchange coefficient for each node in the pipe network. The second option estimates the exchange coefficient using pipe geometry and conduit wall permeability (Equation (13)).
where K j,i,k is the conduit wall permeability term (L T −1 ) in the matrix cell j,i,k, π is the mathematical constant pi (unitless), d ip is the diameter (L) of pipe ip connected to node in, ∆l ip is the straight-line length between nodes (L) of pipe ip connected to node in, τ ip is the tortuosity (unitless) of pipe ip connected to node in, r ip is the radius of the pipe (L). Therefore, α ex is a lumped geometrical and hydraulic property (diameter of the conduit, tortuosity, and the conduit wall permeability) of the karst system [19].

The Geometry of the Conceptual Model
To examine the effect of the exchange flow on the spring hydrograph, a synthetic karst aquifer is proposed. The dimensions and geometry of the synthetic model are analogous to the characteristics of hypothetical models used in previous research [15,19,20,37].
The synthetic model consists of a single conduit with a length of 1700 m which is located in a two-dimensional matrix with a length and width of 1750 and 745 m, respectively (

Types of Recharge
Net recharge (diffuse or concentrated) is considered only from the top of the aquifer. The diffuse recharge is only assigned to the matrix cells. The Gaussian formula is used to generate rainfall events with an interval of 5 min as follows:  The aquifer is considered to be an unconfined layer with a thickness of 100 m. A value of 5 × 10 −3 is assigned to the specific yield. The bottom of the aquifer is considered to be at zero elevation. An initial value of 1 × 10 −5 m/s is assigned to the hydraulic conductivity.
Neumann type (no flow) boundary is defined for the boundaries of the model. The karst conduit is located in the middle of the matrix (10th row) with an elevation of 10 m and discretized in 33 successive pipes with a length of 50 m. A fixed head of 20 m is defined in node 34 at the end of the conduit as the discharge point of the model (Figure 3).
The general flow direction is from the right to the left of the model. The karst conduit is saturated and under pressure in all time steps. The initial hydraulic head is set to 20 m and the water temperature is set at 20°C. The total simulation period of the model is 50,000 min and each step is set to 5 min [37]. The model is initiated in a steady-state at time step 1 (i.e., t = 0).

Types of Recharge
Net recharge (diffuse or concentrated) is considered only from the top of the aquifer. The diffuse recharge is only assigned to the matrix cells. The Gaussian formula is used to generate rainfall events with an interval of 5 min as follows: where f (t) is the cumulative rainfall at each time step (t) (mm), m, δ 2 , a, and b are the constants of 90, 1.5, 20, and 0.12, respectively [37]. m is determinative of the location of the center of the peak, δ 2 together with a determine the extent of kurtosis or flattening, and b is the intensity of rainfall. The rainfall period is assumed from 35 to 210 min from the initial time of the model run.

Discharge Point of the Synthetic Model
Spring is generally considered as the discharge point of the conceptual model. In the designed synthetic model, the spring is directly supplied by the conduit [37,51]. If the spring discharge is mainly supplied by the point recharge (i.e., quick-flow) and the recession curve follows a linear slope (β), the flow regime is dominantly the CFR. However, the flow regime corresponds to the CIFR and/or the MRFR when the spring discharge is mainly supported by the exchange flow from the matrix to the conduit (as the baseflow), consequently, the recession curve follows an exponential slope (α).

Definition of the Scenarios
The direction of the exchange flow between the matrix and the conduit depends on the direction of the hydraulic head difference between the matrix and conduit. Two sets of scenarios including group A, where Q MC+ discharge (the direction of the exchange flow from the matrix to the conduit) was established during all time steps, and group B, where Q MC− discharge (the direction of the exchange flow from the conduit to the matrix) is in early time steps during a rainfall event and Q MC+ discharge is in the next time steps. The characteristics of the scenarios are presented in Table 1. Table 1. The detailed specifications of the scenarios (K is the hydraulic conductivity of the matrix, d is conduit diameter, k c is roughness, S y is specific yield, S S is specific storage, and h initial is the initial hydraulic head in the matrix). (20 m) The first (i.e., A1 and B1), second (i.e., A2 and B2), and third (i.e., A3 and B3) scenarios were designed to examine the effect of the exchange coefficient, the conduit diameter, and the matrix hydraulic conductivity on the shape of the spring hydrograph, respectively.

Number of Run Models Variable Parameter Constant Parameters Type of Scenario
Although the different values of 10 −2 , 10 −3 , 10 −4 , and 10 −5 m 2 /s have been assigned to the exchange coefficient in the A1 and B1 scenarios, the conduit diameter and the matrix hydraulic conductivity remained constant with a value of 0.5 m and 10 −5 m/s, respectively. Moreover, scenarios A2 and B2 were executed by assigning the different conduit diameters (0.1, 0.2, 0.3, 0.4, 0.5 m). The matrix hydraulic conductivity was assumed to be variable (10 −3 , 10 −4 , 10 −5 , 10 −6 m/s) in the third scenarios A3 and B3. In the second and third scenarios, the exchange coefficient was assumed to be constant as 10 −3 m 2 /s.
Based on the direction of exchange flow in scenarios, the flow regime is divided into two groups including (1) CFR, where the direction of exchange flow is from the conduit to the matrix during a restricted part of CFR and it is reversed for the rest of the time and (2) CIFR and/or the MRFR, where the direction of exchange flow is from the matrix to the conduit.
The CRCH Package in CFP is used to define the point recharge (i.e., concentrated recharge) in the two sets of scenarios. The point recharge in nodes 27, 28, 29, and 30, (Figure 3) is equal to the rate of the matrix diffuse recharge in one set (Set A in Table 1) while it is 500 times larger in another set of scenarios (Set B in Table 1).

The Effect of the Exchange Coefficient on the Karst Spring Hydrograph
The results of scenarios A1 and B1 are presented in Figure 4, Supplementary Tables S2 and S3, respectively. Considering scenario A1, the flow regime is mainly controlled by the CIFR and/or the MRFR, as well as by the direction of the exchange flow being from the matrix to the conduit (Q MC+ discharge). Accordingly, the hydrograph recession curve follows an exponential slope ( Figure 4A). However, assuming exchange coefficient value of 10 −5 m 2 /s, the rate of Q MC+ discharge becomes less than the point recharge (Supplementary  Table S2). Therefore, the CFR is established during the earlier time of the recession curve until 130 min. As a result, the slope of the recession curve follows the linear and exponential slopes (both show Q MC+ discharge) for the 90 to 130 min period and 130 min to the end of the simulation time, respectively ( Figure 4A and Supplementary Table S2). As the exchange coefficient increases, the recession curve illustrates segments with steeper slopes (α) as well as a wider range of variations. For example, due to the assigning of 10 −2 m 2 /s to the exchange coefficient, four segments with α of 6.9 × 10 −3 , 1.3 × 10 −3 , 3.0 × 10 −4 , and 5.6 × 10 −5 min −1 were distinguished. However, considering lower values for the exchange coefficient such as 10 −5 m 2 /s, two segments with α values of 7.1 × 10 −3 and 3.6 × 10 −6 min −1 were extracted ( Figure 4A). The results revealed that the decrease in the exchange coefficient from 10 −2 to 10 −5 m 2 /s caused the decrease in the volume of Q MC+ discharge (i.e., the baseflow volume) for the same period from the second inflection point to 50,000 min from 2299.20 to 436.08 m 3 as well as Qp from 1.12 to 0.044 m 3 /min.
Considering scenario B1, the exchange flow from the conduit to the matrix (Q MC− discharge) was only established in part of the CFR time during rainfall in the karst aquifer ( Figure 4B and Supplementary Table S3). The duration of the period corresponds to the time interval where Q MC− discharge depends on the exchange coefficient. Water flux intensity from the conduit to the matrix decreased with the decrease in the exchange coefficient. Moreover, the hydraulic head in the matrix and the conduit are not rapidly balanced with the decrease in the exchange flow. Therefore, the duration of Q MC− discharge becomes longer under CFR conditions. When the exchange coefficient decreased from 10 −2 to 10 −5 m 2 /s, the value of the linear slope of β increased from 0.23 to 0.32 m 3 /min 2 , and the peak discharge increased from 15.14 to 19.49 m 3 /min. However, following scenario A1, the discharge volume of the baseflow for the same period decreased from 2431.00 to 450.34 m 3 and also the number of sub-segments with different slopes, or in other words, the variation range of α (slope of recession curve in the baseflow segment), decreased. In general, the range of variation of α in scenario B1 is higher than the range of variation observed considering scenario A1. In scenario B1, the slope of α was computed to be 2.4 × 10 −2 , 3.9 × 10 −3 , 1.1 × 10 −3 , 3.3 × 10 −4 , and 5.6 × 10 −5 min −1 by assuming exchange coefficient value of 10 −2 m 2 /s ( Figure 4B) as well as 2.3 × 10 −2 and 3.6 × 10 −6 min −1 for exchange coefficient value of 10 −5 m 2 /s ( Figure 4B). Considering both A1 and B1 scenarios, there is no evidence to support that peak time depends on the exchange coefficient. Considering scenario B1, the exchange flow from the conduit to the matrix (QMC− discharge) was only established in part of the CFR time during rainfall in the karst aquifer ( Figure 4B and Supplementary Table S3). The duration of the period corresponds to the time interval where QMC− discharge depends on the exchange coefficient. Water flux intensity from the conduit to the matrix decreased with the decrease in the exchange coefficient. Moreover, the hydraulic head in the matrix and the conduit are not rapidly balanced with the decrease in the exchange flow. Therefore, the duration of QMC− discharge becomes longer under CFR conditions. When the exchange coefficient decreased from 10 −2 to 10 −5 m 2 /s, the value of the linear slope of β increased from 0.23 to 0.32 m 3 /min 2 , and the peak discharge increased from 15.14 to 19.49 m 3 /min. However, following scenario A1, the discharge volume of the baseflow for the same period decreased from 2431.00 to 450.34 m 3 and also the number of sub-segments with different slopes, or in other words, the variation range of α (slope of recession curve in the baseflow segment), decreased. In general, the range of variation of α in scenario B1 is higher than the range of variation observed considering scenario A1. In scenario B1, the slope of α was computed to be 2.4 × 10 −2 , 3.9 × 10 −3 , 1.1 × 10 −3 , 3.3 × 10 −4 , and 5.6 × 10 −5 min −1 by assuming exchange coefficient value of 10 −2 m 2 /s ( Figure 4B) as well as 2.3 × 10 −2 and 3.6 × 10 −6 min −1 for exchange coefficient value of 10 −5 m 2 /s ( Figure 4B). Considering both A1 and B1 scenarios, there is no evidence to support that peak time depends on the exchange coefficient.
Since the water level in the conduits increases due to the heavy point recharge of conduits, QMC− discharge can potentially increase in the early time of the recession curve. The contribution of QMC− discharge in the recession curve is evidenced by the CFR and the Since the water level in the conduits increases due to the heavy point recharge of conduits, Q MC− discharge can potentially increase in the early time of the recession curve. The contribution of Q MC− discharge in the recession curve is evidenced by the CFR and the linear slope in the recession curve. The results revealed that Q MC− discharge is not a long-term process and the direction of exchange flow is reversed in response to the decreasing of the water level in the conduit, although, the flow regime remains as the CFR. After completion of the CFR, the flow is controlled by the CIFR and/or the MRFR and the recession curve follows an exponential slope. Figure 5 is an example of Q MC− discharge that contributes only in the early time of the recession curve in the spring hydrograph assuming the exchange coefficient of 10 −3 m 2 /s. The short-term effect of other phenomena on the spring hydrograph resulting from a temporary increase in the water level in the conduit is investigated by Chang et al. [37]. They showed the turbulent conduit could also strongly influence the early time of the spring hydrograph from the point recharge and may also disappear at the late time of the spring hydrograph.
sion curve follows an exponential slope. Figure 5 is an example of QMC− discharge that contributes only in the early time of the recession curve in the spring hydrograph assuming the exchange coefficient of 10 −3 m 2 /s. The short-term effect of other phenomena on the spring hydrograph resulting from a temporary increase in the water level in the conduit is investigated by Chang et al. [37]. They showed the turbulent conduit could also strongly influence the early time of the spring hydrograph from the point recharge and may also disappear at the late time of the spring hydrograph. Based on the results of scenarios A1 and B1, when one observes the occurrence of a QMC− discharge event during a restricted part of the recession curve, the peak discharge increases when the exchange coefficient decreases. However, the volume of baseflow (from the second inflection point to 50,000 min) decreased. This function can be explained by the fact that during the flow condition where QMC− discharge occurs, a part of the point recharge is transferred to the matrix. This volume of water is drained back to the conduit when the exchange flow goes from the matrix to the conduit (QMC+ discharge). Therefore, with the decrease in the exchange coefficient, the intensity of QMC− discharge decreases, and most of the conduit flow is transferred to the spring as quick-flow.
In contrast, when the flow condition is dominated by the QMC+ discharge over the whole recession curve, the spring discharge is supplied by the mixture of the QMC+ discharge and the point recharge to the conduits. Therefore, spring discharge increases as the exchange coefficients increase. Due to the extra recharge of the matrix under the condition of the QMC− discharge, the volume of the baseflow in scenario B1 is larger than that in scenario A1 where the QMC+ discharge is preponderant over the whole recession curve.

The Influence of the Conduit Diameter on the Karst Spring Hydrograph
The  Based on the results of scenarios A1 and B1, when one observes the occurrence of a Q MC− discharge event during a restricted part of the recession curve, the peak discharge increases when the exchange coefficient decreases. However, the volume of baseflow (from the second inflection point to 50,000 min) decreased. This function can be explained by the fact that during the flow condition where Q MC− discharge occurs, a part of the point recharge is transferred to the matrix. This volume of water is drained back to the conduit when the exchange flow goes from the matrix to the conduit (Q MC+ discharge). Therefore, with the decrease in the exchange coefficient, the intensity of Q MC− discharge decreases, and most of the conduit flow is transferred to the spring as quick-flow.
In contrast, when the flow condition is dominated by the Q MC+ discharge over the whole recession curve, the spring discharge is supplied by the mixture of the Q MC+ discharge and the point recharge to the conduits. Therefore, spring discharge increases as the exchange coefficients increase. Due to the extra recharge of the matrix under the condition of the Q MC− discharge, the volume of the baseflow in scenario B1 is larger than that in scenario A1 where the Q MC+ discharge is preponderant over the whole recession curve. The results of scenario A2 (i.e., the CIFR and/or the MRFR and Q MC+ discharge) revealed that the decrease in the conduit diameter value from 0.5 to 0.1 m caused a decrease in the peak discharge from 0.55 to 0.098 m 3 /min. However, for the same period from the second inflection point to 50,000 min, the baseflow volume increased from 2312.65 to 2325.97 m 3 as the conduit diameter decreased from 0.5 to 0.2 m ( Figure 6A and Supplementary Table S4). Also, these results are observed in the study by Chang et al. [37] in the karst aquifer influenced by the diffuse recharge. Moreover, α also decreased in response to conduit diameter decrease. However, with the change in the conduit diameter from 0.5 to 0.1 m, the time of peak discharge remains constant.
in the peak discharge from 0.55 to 0.098 m 3 /min. However, for the same period from the second inflection point to 50,000 min, the baseflow volume increased from 2312.65 to 2325.97 m 3 as the conduit diameter decreased from 0.5 to 0.2 m ( Figure 6A and Supplementary Table S4). Also, these results are observed in the study by Chang et al. [37] in the karst aquifer influenced by the diffuse recharge. Moreover, α also decreased in response to conduit diameter decrease. However, with the change in the conduit diameter from 0.5 to 0.1 m, the time of peak discharge remains constant. The results showed both QMC− discharge (only in the early time of the recession segment) and QMC+ discharge established in the scenario B2 ( Figure 6B and Supplementary  Table S5). Based on the results, when the conduit diameter decreased from 0.5 to 0.1 m, the peak discharge decreased from 16.90 to 0.62 m 3 /min. However, for the same period from the second inflection point to 50,000 min, the baseflow increased from 2411.30 to 2813.99 m 3 as the conduit diameter decreased from 0.5 to 0.1 m.
The duration of the occurrence of QMC− discharge became longer for the smaller conduit diameters so that for a conduit diameter of 0.1 m, QMC− discharge was dominant for the entire time of the CFR. When conduit diameter decreased from 0.5 to 0.1 m, the slope of the recession curve in the CFR (β) decreased from 0.29 to 0.003 m 3 /min 2 and α also decreased in response to a conduit diameter decrease. The time of the peak discharge increased from 95 to 105 min for the conduit diameters of 0.5 and 0.1 m, respectively.
A smaller conduit diameter can create a larger volume of baseflow. This is due to the larger exchange flow as a result of the more important difference of the hydraulic head The results showed both Q MC− discharge (only in the early time of the recession segment) and Q MC+ discharge established in the scenario B2 ( Figure 6B and Supplementary  Table S5). Based on the results, when the conduit diameter decreased from 0.5 to 0.1 m, the peak discharge decreased from 16.90 to 0.62 m 3 /min. However, for the same period from the second inflection point to 50,000 min, the baseflow increased from 2411.30 to 2813.99 m 3 as the conduit diameter decreased from 0.5 to 0.1 m.
The duration of the occurrence of Q MC− discharge became longer for the smaller conduit diameters so that for a conduit diameter of 0.1 m, Q MC− discharge was dominant for the entire time of the CFR. When conduit diameter decreased from 0.5 to 0.1 m, the slope of the recession curve in the CFR (β) decreased from 0.29 to 0.003 m 3 /min 2 and α also decreased in response to a conduit diameter decrease. The time of the peak discharge increased from 95 to 105 min for the conduit diameters of 0.5 and 0.1 m, respectively.
A smaller conduit diameter can create a larger volume of baseflow. This is due to the larger exchange flow as a result of the more important difference of the hydraulic head between the matrix and conduit (Equation (1)). Also, due to the extra recharge of the matrix under the condition of the Q MC− discharge, the volume of the baseflow in scenario B2 is larger than that in scenario A2 where the Q MC+ discharge is preponderant over the whole recession curve.

The Influence of the Matrix Hydraulic Conductivity on the Karst Spring Hydrograph
The scenarios A3 and B3 consisted of a set of cases with a constant value of the exchange coefficient (10 −3 m 2 /s) and conduit diameter (0.5 m), and different values of the matrix hydraulic conductivity (10 −3 , 10 −4 , 10 −5 , 10 −6 m/s).
The results of scenario A3 are presented in Figure 7A and Supplementary Table S6. Based on the results, when the matrix hydraulic conductivity decreased from 10 −3 to 10 −6 m/s, the peak discharge, the volume of baseflow for the same period from the second inflection point to 50,000 min and the peak discharge time changed from 0.87 to 0.33 m 3 /min, from 2424.35 to 940.250 m 3 , and from 150 to 110 min, respectively. Also, the number of non-linear segments with different slopes of the recession curve increased if a smaller value was assigned to the matrix hydraulic conductivity. As an example, a single α value of 3.5 × 10 −4 and four α values of 7.7 × 10 −3 , 1.3 × 10 −3 , 2.8 × 10 −4 , and 1.3 × 10 −5 min −1 were computed by assigning 10 −3 and 10 −6 m/s to the matrix hydraulic conductivity, respectively.

Generalization of the Effect of Exchange Flow on the Spring Hydrograph
The results of the scenarios (Section 4) are used to develop a conceptual mode behavior of the exchange flow in the karst spring hydrograph (Figure 8). To exp systematic variation of the spring hydrograph in response to the variation of the ex flow, two types of recharge are assumed in Figure 8: the diffuse and point recha infiltrates slowly and quickly into the matrix and the conduit, respectively. The dif in type and rate of the recharge causes QMC− and/or QMC+ discharge due to the dif in the hydraulic head between the matrix and conduit. The change of QMC+ to Q charge is influenced by the flow regime (i.e., CIFR and/or the MRFR and CFR) rate of point recharge. In this conceptual model (Figure 8), the rate of point r changes from low to high while the diffuse recharge is assumed to be at a consta Moreover, the recharge event takes place during the early time and the spring di starts to increase in response to the recharge event.
Based on the history of recharge and exchange flow, it can be expected that ferent values of exchange coefficient, conduit diameter, and matrix hydraulic cond play different roles in creating hydrographs with different characteristics such The results of scenario B3 (including both flow conditions: Q MC− and Q MC+ ) are presented in Figure 7B and Supplementary Table S7. The Q MC− discharge was only established in the early time of the CFR (i.e., the beginning segments of the recession curve). Based on the results, when the matrix hydraulic conductivity decreased from 10 −3 to 10 −6 m/s, the duration of occurrence of Q MC− discharge was estimated shorter, the peak discharge increased from 16.27 to 17.54 m 3 /min, and β increased from 0.27 to 0.29 m 3 /min 2 . However, the volume of baseflow for the same period from the second inflection point to 50,000 min decreased from 2600.20 to 980.55 m 3 . The rate of the exponential slope (α) in the segments with the CIFR and/or the MRFR increased in response to assigning a smaller value to matrix hydraulic conductivity. For example, the exponential slope (α) values of 2.6 × 10 −2 , 4.8 × 10 −4 , and 3.5 × 10 −4 min −1 and 5.3 × 10 −2 , 3.9 × 10 −3 , 1.2 × 10 −3 , 2.8 × 10 −4 , and 1.2 × 10 −5 min −1 were estimated when the matrix hydraulic conductivity considered to be 10 −3 and 10 −6 m/s, respectively. In scenario B3, the time of the peak discharge time remained constant when the matrix hydraulic conductivity changed.
Based on the results of scenarios A3 and B3, the movement of the water entering the matrix due to Q MC− discharge is slow if the small values are assigned to the matrix hydraulic conductivity. As a result, the water head is locally increased in short distances around the conduit. Under such a condition (i.e., assigning a small value to the matrix hydraulic conductivity), the larger volume of the point recharge to the conduit is quickly transferred to the spring and causes the volume of the quick flow and β to increase. On the other hand, the volume of the baseflow decreases for the same period from the second inflection point to 50,000 min.

Generalization of the Effect of Exchange Flow on the Spring Hydrograph
The results of the scenarios (Section 4) are used to develop a conceptual model of the behavior of the exchange flow in the karst spring hydrograph (Figure 8). To explain the systematic variation of the spring hydrograph in response to the variation of the exchange flow, two types of recharge are assumed in Figure 8: the diffuse and point recharge that infiltrates slowly and quickly into the matrix and the conduit, respectively. The difference in type and rate of the recharge causes Q MC− and/or Q MC+ discharge due to the difference in the hydraulic head between the matrix and conduit. The change of Q MC+ to Q MC− discharge is influenced by the flow regime (i.e., CIFR and/or the MRFR and CFR) and the rate of point recharge. In this conceptual model (Figure 8), the rate of point recharge changes from low to high while the diffuse recharge is assumed to be at a constant rate. Moreover, the recharge event takes place during the early time and the spring discharge starts to increase in response to the recharge event.
Based on the history of recharge and exchange flow, it can be expected that the different values of exchange coefficient, conduit diameter, and matrix hydraulic conductivity play different roles in creating hydrographs with different characteristics such as peak discharge and baseflow volume (Figure 8). In this comparison, with the change of each of the parameters of interest (exchange coefficient, conduit diameter, and matrix hydraulic conductivity), other parameters of the karst aquifer are assumed constant.
In column A, there is no CFR due to small amounts of point recharge in the karst aquifer. Under such conditions, the spring hydrograph is influenced by CIFR and/or the MRFR and the direction of the exchange flow is from the matrix to the conduit (Q MC+ discharge). In the CIFR and/or the MRFR, the recession curve follows an exponential (α) equation. Under these conditions, the simulated hydrographs present greater peak discharge assigning a large value to each of the parameters of interest (the exchange coefficient, the conduit diameter, and the hydraulic conductivity of the matrix). The number of segments of the recession curve with different slopes increase when the exchange coefficient is large. However, they decrease when the hydraulic conductivity of the matrix is large. In column A, there is no CFR due to small amounts of point recharge in the karst aquifer. Under such conditions, the spring hydrograph is influenced by CIFR and/or the Also, the inflection point in the recession curve indicates the end of the infiltration and recharge process. In column B, although the CFR can be established in the early time steps of hydrograph by increasing the amount of point recharge, exchange flow direction may still be from the matrix to the conduit (Q MC+ discharge). In the CFR, the recession curve follows a linear relationship with the β-slope. Under these conditions, the inflection point in the recession curve indicates the end of the CFR as well as the end of the infiltration and recharge processes. Similar to column A, under these conditions, the higher values of each of the parameters of interest (the exchange coefficient, the conduit diameter, and the hydraulic conductivity of the matrix) make the larger peak discharge. Finally, in column C, the direction of the exchange flow changes from the conduit to the matrix (Q MC− discharge) when the point recharge increases significantly. In this case, the spring hydrograph experiences three modes including (1) CFR and Q MC− discharge in the early times of the hydrograph, (2) the CFR and Q MC+ discharge, and (3) the CIFR and/or the MRFR and Q MC+ discharge in the rest of the time. In modes of 1 and 2, the recession curve follows a linear relationship with the β-slope, and in the mode of 3, the recession curve follows an exponential relationship with the α-slope. Then, the Q MC− discharge and conduit flow regime consist of short-term processes. Over time, the flow regime changes to CIFR and/or the MRFR, as well as the exchange flow direction, which changes from the matrix to the conduit because the drop of the hydraulic head in the conduit is faster than that in the matrix. Under the conditions indicated in column C, where the karst system experiences Q MC− discharge in the part of CFR due to a large point recharge, the response of the spring hydrograph to the exchange coefficient, the conduit diameter, and the hydraulic conductivity of the matrix is different from that in columns A and B. Therefore, the simulated hydrographs present greater peak discharge when the exchange coefficient is small. As a result, the small volume of water is stored in the matrix due to the exchange flow and the majority of water enters into a conduit network that is drained by the spring. Moreover, under the condition that constitutes Q MC− discharge (column C), the volume of baseflow will be larger than the condition when the karst system has only experienced Q MC+ discharge (columns A and B) for a specific exchange coefficient. This phenomenon is evidenced by a flow pipe as spring outflow volume (the higher the volume of the discharge caused by the larger value of the pipe diameter). Based on the results of a conceptual modeling approach, the influence of matrix hydraulic conductivity on the spring hydrograph is the same as the role of the exchange coefficient in different flow regimes and directions of exchange flow. However, regardless of the exchange flow direction (from the matrix to the conduit or conversely) that occurs in the karst aquifer, larger conduit diameters produce hydrographs with larger discharge peaks. However, a smaller conduit diameter can create a larger volume of baseflow. This is due to the larger exchange flow as a result of the more important difference in the hydraulic head between the matrix and conduit.
To classify spring hydrographs and conclude for possible internal characteristics of a karst aquifer, a systematic approach is proposed based on the reliable concluding points obtained from the analysis of the shape of the synthetic spring hydrographs that resulted from scenarios A and B (Table 2). Based on the defined scenarios, the amount of the conduit flow in scenario B is more important than that in scenario A, although, the rate of recharge into the matrix is the same in both scenarios. The matrix properties such as the specific yield and the aquifer thickness are constant in both scenarios for the aquifer with the no-flow boundaries where the outflow of the matrix is through the exchange flow from the matrix to the conduit and then spring. the no-flow boundaries where the outflow of the matrix is through the exchange flow from the matrix to the conduit and then spring. The effective porosity of the matrix in Scenario A = The effective porosity of the matrix in Scenario B αex: Exchange coefficient T: Transmissivity α: The slope of recession curve under the condition of baseflow Δh: The difference of the hydraulic head between the matrix and the conduit RCH: Diffuse recharge to the matrix CRCH: Concentrated recharge to conduits Given a spring hydrograph, following the proposed approach in Table 2 helps to assess preliminary knowledge about the exchange flow as well as the internal karst structure. Firstly, the logarithmic recession curve should be derived based on the original spring hydrograph ( Table 2). The visual and clear existence or lack of inflection point on the recession curve is a decision point. If the inflection point is visually clear, its position should be identified. Under this condition, 3 modes can be imagined.
In mode І, the inflection point is located before an arc region. This condition is observed in scenarios A and B by assigning large values to the exchange coefficient and small values to the transmissivity of the matrix.
In modes Ⅱ and Ⅲ , the inflection point is located before a non-arc region that represents a rapid transition from the flood recession to the baseflow recession. Mode Ⅱ is  The effective porosity of the matrix in Scenario A = The effective porosity of the matrix in Scenario B αex: Exchange coefficient T: Transmissivity α: The slope of recession curve under the condition of baseflow Δh: The difference of the hydraulic head between the matrix and the conduit RCH: Diffuse recharge to the matrix CRCH: Concentrated recharge to conduits Given a spring hydrograph, following the proposed approach in Table 2 helps to assess preliminary knowledge about the exchange flow as well as the internal karst structure. Firstly, the logarithmic recession curve should be derived based on the original spring hydrograph ( Table 2). The visual and clear existence or lack of inflection point on the recession curve is a decision point. If the inflection point is visually clear, its position should be identified. Under this condition, 3 modes can be imagined.
In mode І, the inflection point is located before an arc region. This condition is observed in scenarios A and B by assigning large values to the exchange coefficient and small values to the transmissivity of the matrix.
In modes Ⅱ and Ⅲ , the inflection point is located before a non-arc region that represents a rapid transition from the flood recession to the baseflow recession. Mode Ⅱ is  The effective porosity of the matrix in Scenario A = The effective porosity of the matrix in Scenario B αex: Exchange coefficient T: Transmissivity α: The slope of recession curve under the condition of baseflow Δh: The difference of the hydraulic head between the matrix and the conduit RCH: Diffuse recharge to the matrix CRCH: Concentrated recharge to conduits Given a spring hydrograph, following the proposed approach in Table 2 helps to assess preliminary knowledge about the exchange flow as well as the internal karst structure. Firstly, the logarithmic recession curve should be derived based on the original spring hydrograph ( Table 2). The visual and clear existence or lack of inflection point on the recession curve is a decision point. If the inflection point is visually clear, its position should be identified. Under this condition, 3 modes can be imagined.
In mode І, the inflection point is located before an arc region. This condition is observed in scenarios A and B by assigning large values to the exchange coefficient and small values to the transmissivity of the matrix.
In modes Ⅱ and Ⅲ , the inflection point is located before a non-arc region that represents a rapid transition from the flood recession to the baseflow recession. Mode Ⅱ is  The effective porosity of the matrix in Scenario A = The effective porosity of the matrix in Scenario B αex: Exchange coefficient T: Transmissivity α: The slope of recession curve under the condition of baseflow Δh: The difference of the hydraulic head between the matrix and the conduit RCH: Diffuse recharge to the matrix CRCH: Concentrated recharge to conduits Given a spring hydrograph, following the proposed approach in Table 2 helps to assess preliminary knowledge about the exchange flow as well as the internal karst structure. Firstly, the logarithmic recession curve should be derived based on the original spring hydrograph ( Table 2). The visual and clear existence or lack of inflection point on the recession curve is a decision point. If the inflection point is visually clear, its position should be identified. Under this condition, 3 modes can be imagined.
In mode І, the inflection point is located before an arc region. This condition is observed in scenarios A and B by assigning large values to the exchange coefficient and small values to the transmissivity of the matrix.
In modes Ⅱ and Ⅲ , the inflection point is located before a non-arc region that represents a rapid transition from the flood recession to the baseflow recession. Mode Ⅱ is  The effective porosity of the matrix in Scenario A = The effective porosity of the matrix in Scenario B αex: Exchange coefficient T: Transmissivity α: The slope of recession curve under the condition of baseflow Δh: The difference of the hydraulic head between the matrix and the conduit RCH: Diffuse recharge to the matrix CRCH: Concentrated recharge to conduits Given a spring hydrograph, following the proposed approach in Table 2 helps to assess preliminary knowledge about the exchange flow as well as the internal karst structure. Firstly, the logarithmic recession curve should be derived based on the original spring hydrograph ( Table 2). The visual and clear existence or lack of inflection point on the recession curve is a decision point. If the inflection point is visually clear, its position should be identified. Under this condition, 3 modes can be imagined.
In mode І, the inflection point is located before an arc region. This condition is observed in scenarios A and B by assigning large values to the exchange coefficient and small values to the transmissivity of the matrix.
In modes Ⅱ and Ⅲ , the inflection point is located before a non-arc region that represents a rapid transition from the flood recession to the baseflow recession. Mode Ⅱ is  Given a spring hydrograph, following the proposed approach in Table 2 helps to assess preliminary knowledge about the exchange flow as well as the internal karst structure. Firstly, the logarithmic recession curve should be derived based on the original spring hydrograph ( Table 2). The visual and clear existence or lack of inflection point on the recession curve is a decision point. If the inflection point is visually clear, its position should be identified. Under this condition, 3 modes can be imagined.
In mode І, the inflection point is located before an arc region. This condition is observed in scenarios A and B by assigning large values to the exchange coefficient and small values to the transmissivity of the matrix.
In modes Ⅱ and Ⅲ , the inflection point is located before a non-arc region that represents a rapid transition from the flood recession to the baseflow recession. Mode Ⅱ is the no-flow boundaries where the outflow of the matrix is through the exchange flow from the matrix to the conduit and then spring. Given a spring hydrograph, following the proposed approach in Table 2 helps to assess preliminary knowledge about the exchange flow as well as the internal karst structure. Firstly, the logarithmic recession curve should be derived based on the original spring hydrograph ( Table 2). The visual and clear existence or lack of inflection point on the recession curve is a decision point. If the inflection point is visually clear, its position should be identified. Under this condition, 3 modes can be imagined.
In mode І, the inflection point is located before an arc region. This condition is observed in scenarios A and B by assigning large values to the exchange coefficient and small values to the transmissivity of the matrix.
In modes Ⅱ and Ⅲ , the inflection point is located before a non-arc region that represents a rapid transition from the flood recession to the baseflow recession. Mode Ⅱ is  Given a spring hydrograph, following the proposed approach in Table 2 helps to assess preliminary knowledge about the exchange flow as well as the internal karst structure. Firstly, the logarithmic recession curve should be derived based on the original spring hydrograph ( Table 2). The visual and clear existence or lack of inflection point on the recession curve is a decision point. If the inflection point is visually clear, its position should be identified. Under this condition, 3 modes can be imagined.
In mode І, the inflection point is located before an arc region. This condition is observed in scenarios A and B by assigning large values to the exchange coefficient and small values to the transmissivity of the matrix.
In modes Ⅱ and Ⅲ , the inflection point is located before a non-arc region that represents a rapid transition from the flood recession to the baseflow recession. Mode Ⅱ is Given a spring hydrograph, following the proposed approach in Table 2 helps to assess preliminary knowledge about the exchange flow as well as the internal karst structure. Firstly, the logarithmic recession curve should be derived based on the original spring hydrograph ( Table 2). The visual and clear existence or lack of inflection point on the recession curve is a decision point. If the inflection point is visually clear, its position should be identified. Under this condition, 3 modes can be imagined.
In mode I, the inflection point is located before an arc region. This condition is observed in scenarios A and B by assigning large values to the exchange coefficient and small values to the transmissivity of the matrix.
In modes I and III, the inflection point is located before a non-arc region that represents a rapid transition from the flood recession to the baseflow recession. Mode II is observed in scenarios A and B by assigning small values to the exchange coefficient and different large and small values to the transmissivity of the matrix. Under these conditions, the conduit flow is more important than that in mode I.
The inflection point is not distinguishable in model IV (Table 2). By assigning large values to the exchange coefficient and the transmissivity of the matrix, the obtained hydrographs were similar to mode III in the scenario B and mode IV in scenario A and the inflection point on the logarithmic curve of the recession was not observed (Table 2).
A real field example of this condition is found in the Santa Fe River Rise, studied by Bailly-Comte et al. [16].
Based on the result of the modeling for the synthetic karst aquifer, when the exchange coefficient is large enough, the slope of the recession curve under the condition of the baseflow is higher than when the exchange coefficient is small (the modes of I, III, and IV). Also, by assigning small values to the exchange coefficient (mode II), the difference of the hydraulic head between the conduit and the matrix is more important than when large values are assigned to the exchange coefficient.
Despite extensive attempts in this research including the running of a lot of simulations under the different characteristics of the synthetic karst aquifer such as the matrix with different dimensions, saturation thickness, different rate of recharge, and effective porosity, it is unlikely to introduce a global value for the large/small values of the exchange coefficient and the difference of the hydraulic head between the conduit and the matrix in karst aquifers. It seems that the large/small values of the exchange flow depend on the conditions of each karst aquifer. Therefore, a definite value cannot be proposed as the large/small values for all aquifers. This topic requires more detailed research and can be suggested as an open question for further research.

Conclusions
The main focus of the study was to investigate the effect of the matrix-conduit exchange flow (amount and direction) on the karst spring hydrograph. Assuming different values of the exchange coefficient and the hydraulic head difference between the matrix and the conduit allow for the investigation of the effect of matrix-conduit exchange flow magnitude on a spring hydrograph. The difference of hydraulic head between the matrix and the conduit influences the direction of exchange that in turn depends on the amount and type of recharge (point recharge or diffuse recharge). In most cases, fluxes go from the matrix to the conduit. Increasing the point recharge leads to a gradient inversion between matrix and conduit, as well as the directions of the exchange flow. After a rainfall event, the water entering the matrix under a short-term process of the Q MC− discharge returns to the conduit as baseflow. Under this condition, depending on the gradient and the magnitude of the exchange flow, the volume of baseflow increases and becomes greater than the scenario where there is only Q MC+ discharge during recharge and recession events (as seen in scenario A). The shape of the spring hydrograph can be influenced by the direction of the exchange flow during the recharge and recession events. Hence, if the direction of exchange flow is from the conduit to the matrix, even during a restricted part of the conduit flow regime, the simulated hydrographs present greater peak discharge when each of the parameters of interest (the exchange coefficient or the hydraulic conductivity of the matrix) are small.
According to the different scenarios, it appears that the amount of exchange flow may strongly influence the spring discharge of the karst aquifers accordingly in the karst system conceptualization. However, there is a challenge to estimate the matrix-conduit interaction since it is difficult and quite impossible to quantify in the field.
Assessing the matrix-conduit dynamics can be used for a wide range of applications in real karst systems. The following issues can be investigated concerning the matrixconduit dynamics: Is aquifer vulnerability to contaminants influenced by matrix-conduit dynamics based on both the short-term and long-term?
Can a punctual infiltration of contaminant deteriorate the capacitive function of an aquifer?
May the groundwater abstraction influence the dynamics of the quick and slow flow? Hence, deteriorating the local water resource? In particular, in the coastal aquifers under pumping, it can be one of the major environmental concerns. Geng and Michael [52] and Kreyns et al. [53] showed the importance of conduits and preferential flow paths on the complexity of the saltwater intrusion distribution pattern and mixing dynamics between saline and fresh groundwater.
This can also lead to better constrain the karst structure (i.e., conduit geometry). This kind of research may also be considered together with solute transport modeling. Therefore, the correct estimation of the exchange flow (amount and direction) is of prior importance when dealing with the above-mentioned issues in a real case study.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/w13091189/s1: Table S1. Reviews of list of introduced formulas for computing the exchange flow. Table S2. The characteristics of the spring hydrograph in the scenario A1. Table S3. The characteristics of the spring hydrograph in scenario B1. Table S4. The characteristics of spring hydrograph in response to changes of conduit diameter in scenario A2. Table S5. The characteristics of spring hydrograph in response to changes in conduit diameter in scenario B2. Table S6. The spring hydrograph changes influenced by matrix hydraulic conductivity changes in scenario A3. Table S7. The spring hydrograph changes influenced by matrix hydraulic conductivity changes in scenario B3.