Representation of water abstraction from a karst conduit with numerical discrete-continuum models

. Karst aquifers are characterized by highly conductive conduit ﬂow paths embedded in a less conductive ﬁssured and fractured matrix, resulting in strong permeability contrasts with structured heterogeneity and anisotropy. Groundwater storage occurs predominantly in the ﬁs-sured matrix. Hence, most mathematical karst models assume quasi-steady-state ﬂow in conduits neglecting conduit-associated drainable storage (CADS). The concept of CADS considers storage volumes, where karst water is not part of the active ﬂow system but hydraulically connected to conduits (for example karstic voids and large fractures). The disregard of conduit storage can be inappropriate when direct water abstraction from karst conduits occurs, e.g., large-scale pumping. In such cases, CADS may be relevant. Furthermore, the typical ﬁxed-head boundary condition at the karst outlet can be inadequate for water abstraction scenarios because unhampered water inﬂow


Introduction
Karst aquifers can be described as triple porosity systems with continuous primary porosity in the matrix, secondary porosity within fissures and fractures, and tertiary porosity represented by solution-enlarged features, i.e., highly permeable conduits (Worthington et al., 2000). Different conceptual approaches regarding karst aquifer storage are presented in literature (Bakalowicz, 2005). Mangin (1975Mangin ( , 1994 attributed large water storage to poorly interconnected large voids adjacent to conduit systems, which transmit water from the groundwater table towards a spring (annex systems to drainage). Mangin (1994) further assumes that matrix storage is negligible. In contrast, Drogue (1974Drogue ( , 1992 proposed areal water storage in the hydraulically continuous matrix drained by the highly permeable karst conduit system. Storage directly associated with discrete conduits was not considered. However, the existence of water storage within highly permeable karst structures is known from karst hydraulics (e.g., Cornaton and Perrochet, 2002;Geyer et al., 2008;Maréchal et al., 2008). Worthington et al. (2000) and Worthington (2007) presented field studies that clearly emphasize the necessity to describe karst aquifers as triple porosity systems.

Published by Copernicus Publications on behalf of the European Geosciences Union.
Strong permeability contrasts within triple porosity systems lead to uncertainties during characterization of these systems with conventional experimental, analytical, and numerical methods (Geyer et al., 2013). For example, artificial tracer tests are suitable methods for characterization of large karst conduit systems (e.g., Atkinson et al., 1973). However, tracer methods fail to characterize matrix properties on a catchment scale. In contrast, conventional hydraulic borehole tests provide useful information about local aquifer properties but also are unable to determine hydraulic information on catchment scale because of limited pumping rates.
Only a few experiments are documented that address an integral large-scale characterization of karst aquifers. In particular, Maréchal et al. (2008) performed a pumping test with abstraction rates up to several hundred liters per second for about one month. These high abstraction rates were possible because the pumping well was directly connected to the conduit system. The pumping test produced drawdowns in both the conduit system and the fissured matrix and, therefore, provided evidence to infer the existence of large water storages within the fissured matrix and the conduits (Maréchal et al., 2008), see Fig. 1. Amongst others, diagnostic plots of drawdown and the logarithmic derivative of the drawdown over time on a log-log plot (Bourdet et al., 1983) were used to interpret the pumping test drawdown. At early times, the drawdown and the derivative follow a straight line of unit slope indicating the dominance of storage effects (e.g., Bourdet et al., 1983;Renard et al., 2009), i.e., wellbore storage, or interchangeable for directly pumping from karst conduits, karst conduit storage. This unit slope of drawdown and derivative is present during about the first 1000 min of the large-scale pumping test, i.e., for early times (Fig. 6 in Maréchal et al., 2008).
Further evaluation of large-scale field experiments for karst characterization can be achieved by numerical modeling, ideally with an approach that considers both matrix and conduits with distributed parameter fields. Lumpedparameter modeling approaches for simulation of karst hydraulics do not consider a distributed hydraulic parameter field, e.g., Geyer et al. (2008), Maréchal et al. (2008) and others. Halihan and Wicks (1998) introduced a model approach of pipes with smoothly turbulent flow that connect reservoirs to represent conduit-flow aquifers, i.e., flow into or out of the matrix is not considered. This idea to align reservoirs and flow restrictions in series was applied by Prelovšek et al. (2008) to karst systems in Slovenia. Covington et al. (2009) presented physically more enhanced representations of single karst network elements like full pipes, open channels and reservoirs. However, a primary limitation of this model is that it is only applied to single karst network elements without the consideration of matrix interaction. Discrete-continuum (or "hybrid") models, which couple a discrete pipe flow model to a continuum model, represent a suitable approach to simulate karst aquifers (e.g., Király, 1998Király, , 2002Sauter et al., 2006;Kaufmann, 2009) without neglecting the matrix interaction and under consideration of strongly anisotropic hydraulic parameter fields.  presented a discrete-continuum model for the simulation of laminar and turbulent pipe flow coupled to a continuum that represents the matrix. This approach was further developed as Conduit Flow Process Mode 1 (CFPM1) for MODFLOW-2005(Shoemaker et al., 2008. This discretecontinuum model approach was applied in a number of modeling studies presented in scientific literature (e.g., Liedl et al., 2003;Bauer et al., 2003;Hill et al., 2010). However, discrete-continuum models based on the Liedl et al. (2003) approach simulate conduit flow as quasi-steady without drainable storage. Consequently, these discrete-continuum models fail to simulate transient water storage within the conduit system because steady-state pipe flow equations are applied, i.e., drainable conduit storage is not considered. For this reason, transient water storage can be considered by the continuum model only. Following this, Reimann et al. (2011) presented a discrete-continuum modeling approach to simulate unsteady discrete flow in a variably filled pipe network employing the Saint-Venant equations, while de Rooij (2008) and de Rooij et al. (2013) additionally considered surface flow and variably saturated matrix flow. The approach appears to be suitable for fundamental studies but high parameter demand and computational effort can be identified as potential drawbacks.
The objective of this work is to provide a distributive process-based modeling approach based on the  approach that allows the simulation of hydraulic impacts (water abstraction, large-scale hydraulic tests, discharge events) on karst systems with manageable data demand and computational effort, under consideration of important storage processes. For that reason, the discretecontinuum modeling approach of CFPM1 was further enhanced by adding water storage in parallel to the conduits. In the following, we refer to this as conduit-associated drainable storage (CADS). Scenarios with direct water abstraction from the conduits can result in a reversion of the flow direction, i.e., inflow at the karst spring (e.g., Maréchal et al., 2008). Hence, the numerical model was complemented by a constrained boundary condition according to Bauer et al. (2005) to avoid unhampered water inflow through the spring. The performance of the numerical model is evaluated by a verification test and highly simplified synthetic scenarios. Subsequently, an idealized model representation of the large-scale pumping test at the Cent Fonts karst system" behind large-scale pumping test (Maréchal et al., 2008) is considered as application outlook to demonstrate the potential and benefits of the enhanced discrete-continuum model under field conditions.  (Maréchal et al., 2008); right: abstraction rate and drawdown in both matrix and conduit for a long-term and large-scale pumping test; figures from Maréchal et al. (2008).

Modeling approach
This section introduces the concept of conduit-associated drainable storage and implementation of this concept within the numerical discrete-continuum model CFPM1 (Shoemaker et al., 2008). Further explanation is given about the implementation of a flow constrained boundary condition in CFPM1.

Conceptual consideration of conduit-associated drainable storage
In general, storage in karst systems occurs in (A) the porous matrix (primary porosity), (B) fractures/fissures (secondary porosity), and (C) solution-enlarged pathways like conduits (tertiary porosity), Fig. 2 left. The discrete-continuum model concept considers two compartments: (1) a representative elementary volume (REV) of the fissured/fractured matrix simulated as continuum with laminar flow and storage (matrix continuum), and (2) solution-enlarged highly permeable conduits simulated as discrete pipe network with laminar and turbulent flow without storage (active flow system). Existing discrete-continuum models with steady pipe flow equations provide drainable storage only by the matrix continuum. Dynamic processes like water abstraction, however, demonstrate that additional fast-responding local storage is present (e.g., Maréchal et al., 2008). This fast responding storage is assumed to be provided by features like solutionenlarged fractures (SF1), and other cavities (SF2) that are directly associated (connected) to the conduit flow system but do not actively participate in pipe flow, i.e., conduitassociated drainable storage. Figure 2 illustrates this concept.

Technical implementation into MODFLOW-2005
Conduit Flow Process (CFP) 2.2.1 Numerical discrete-continuum model CFP Mode 1 Shoemaker et al. (2008) incorporated the discrete-continuum model approach of Liedl et al. (2003) in MODFLOW-2005 as CFPM1. Laminar groundwater flow in the continuum model is represented by the Darcy equation:  ald and Harbaugh, 1988). Further details regarding the solution of the groundwater flow equation are well documented in MODFLOW manuals and, therefore, not included here.
The conduit system is represented by nodes that are connected by cylindrical pipes. Detailed explanation about the discrete pipe model is given by Shoemaker et al. (2008). Following, a short overview is provided. Volume conservation at each node is expressed by Kirchhoff's law (Horlacher and  Lüdecke, 1992): where Q ip is discharge from pipe i (up to n pipes) [L 3 T −1 ] and Q ss is the sum of flow from sinks and sources [L 3 T −1 ]. Laminar pipe flow is represented by the Hagen-Poiseuille equation: with d ip diameter of pipe i [L], g gravitational acceleration [LT −2 ], h c,ip head difference along pipe i, ν kinematic viscosity of water [L 2 T −1 ], and l ip length of pipe i [L] (Shoemaker et al., 2008). Turbulent pipe flow is considered by the Darcy-Weisbach equation with the friction factor according to the Colebrook-White equation (Shoemaker et al., 2008): with k c,ip mean roughness height of pipe i [L]. The coupling between pipe network and continuum model is realized through a head-dependent exchange flow rate Q ex : where h c [L] is the conduit head, h m [L] is the matrix head, and α ex [L 2 T −1 ] is the pipe conductance (Barenblatt, 1960;Shoemaker et al., 2008). Water transfer Q ex between matrix and conduits is considered for each node by MODFLOW as external flow and by the pipe system as sink or source (term Q ss in Eq. 2).

Implementation of CADS
To consider drainable conduit storage within CFPM1, the CADS package was developed. Conceptually, CAD storage is assumed to be in direct hydraulic contact with draining conduits, see Fig. 3, so that where h CADS is the head [L] in the CAD storage. The CAD storage is conceptualized as a rectangular block directly associated with the pipe, i.e., the CADS base area is the pipe length associated with a node times the width of the CAD storage (Fig. 3). Finally, the volume of the CADS for each node is computed as where l CADS is the length [L] (= length of pipes associated to a node), W CADS is the width of the CAD storage [L], and z bot is the elevation of the pipe bottom [L]. The ratio V CADS /(h CADSz bot ) corresponds to the free-surface area of the dewatering conduit network defined by Maréchal et al. (2008). It is assumed that water released from the CADS due to head variations immediately enters the pipe resulting in additional discharge. The resulting discharge from CADS storage, Q CADS [L 3 T −1 ], is considered as where V t is the volume of the CAD storage [L 3 ] at the time t and t is the time step size [T]. Q CADS is directly added to the CFPM1 system of equations (represented by Q ss in Eq. 2) and subsequently considered for the iterative solution.

Implementation of a constrained fixed-head boundary condition
A conduit with a fixed-head boundary condition can strongly affect in-or outflow of the highly permeable pipe network at the outlet, i.e., a karst spring. For example, water abstraction from a distinctive and well developed pipe network can result in unlimited water inflow through a fixed-head boundary. However, this contradicts the drawdown behavior in field situations (e.g., Maréchal et al., 2008;Fig. 1). The fixed-head limited-flow (FHLQ) boundary condition is intended to limit inflow for constant head boundaries. If a user-defined discharge threshold is exceeded, the fixed-head boundary condition switches to a fixed-flow boundary condition, which results in a variable head (Bauer et al., 2005): with H fixed head value (FH) [L], Q discharge at the boundary (negative values denote outflow) [L 3 T −1 ], and Q L limiting discharge (LQ) [L 3 T −1 ]. H and Q L are to be defined by the user according to site-specific conditions.

Test scenarios
The functionality of conduit-associated drainable storage in CFPM1 is verified by draining an isolated conduit with CADS. Subsequently, a highly simplified model is used to test the interaction of a conduit with CADS and the matrix continuum for water abstraction setups under specific boundary conditions.

Verification test with an isolated conduit
The test setup considers an l p = 500 m-long pipe that is subdivided by 6 nodes and 5 equally long tubes with a radius of r = 0.05 m and a bottom elevation of 0 m. CAD storage is considered for the upstream node only (node 1) with W CADS = 0.1 m and a node-associated conduit length of l CADS = 50 m. An initial inflow of Q 0 = 1.0 m 3 s −1 is applied to node 1 and a fixed head of 50 m is considered as downstream boundary condition at node 6. Inflow stops immediately at t = 0 and, subsequently, CAD storage is drained. The resulting (drainage) flow Q d can be described by the recession function from Maillet (1905).
with K c hydraulic conductivity of the conduit [LT −1 ]. Conduit flow is assumed as laminar with K c = 2340 ms −1 according to the Hagen-Poiseuille equation (e.g., Shoemaker et al., 2008, p. 8). The resulting upstream head in node 1 is 77.163 m, respectively h is 27.163 m. The recession discharge along time, computed by Eq. (10) and by CFPM1 with CADS, is presented in Fig. 4. Both results are equal and, therefore, demonstrate the ability of CADS to represent the dynamic behavior of storage that is directly coupled with a conduit. After drainage is completed, i.e., conduit heads are equal to the fixed head of 50 m, CFP budget files account for 135.814 m 3 of water released from CADS. This equals W CADS × l CADS × h = 0.1 m × 50 m × 27.163 m and, therefore, verifies the implementation of CADS in CFPM1.

Simple coupled system
In this section, the interactions of pipes with CADS with a matrix continuum under different boundary conditions are investigated. The intention of the test examples is to demonstrate the functioning of the model enhancements in a simplified (and therefore traceable) environment to allow a systematic process study.

Basic model setup
The basic model setup consists of a continuum model with 11 columns and 11 rows where each cell is 100 m × 100 m (Fig. 5). Hydraulic conductivity and storage coefficient of the matrix continuum are set as K m = 1 × 10 −5 ms −1 and S m = 0.01, respectively. The embedded conduit consists of 6 nodes connected by 5 tubes (each 100 m long). The pipe diameter for one model realization is 0.5 m and for another model realization 2.5 m. Pipe roughness height k c is set to 0.01 m. Water transfer between conduits and matrix is parameterized by a fixed water transfer coefficient per unit length α ex /l p = 1 × 10 −5 ms −1 . All lateral outer boundaries of the matrix continuum are of Neumann type (no flow). Diffuse areal recharge is uniformly applied to the continuum with 8.26 × 10 −8 ms −1 and direct point recharge of 0.1 m 3 s −1 is applied to conduit node 1 (Fig. 4). The karst spring is represented by a fixed head of 50 m at conduit node 6. Water abstraction occurs in node 5. The model simulates three periods, specifically: (1) pre-pumping from 0 to 86 400 s (1 day duration); (2) pumping from 86 400 to 345 600 s (3 days duration) at rates of 0.3, 0.5, and 1.0 m 3 s −1 , respectively, and (3) recovery from 345 600 to 604 800 s (3 days duration).

Results for the basic model (available CFPM1 without CADS)
The resulting conduit heads, flow from matrix transfer and flow at the fixed head (node 6) along time are shown in Fig. 6. For stress period 1, spring discharge equaled 0.2 m 3 s −1 (denoted as negative flow at the fixed head) and consisted of diffuse areal recharge (0.1 m 3 s −1 ) that enters the conduit as matrix transfer flow plus direct point recharge to the conduit (0.1 m 3 s −1 at node 1, see Fig. 5). Water abstraction in period 2 from node 5 results in an immediate conduit head drawdown to quasi-steady conduit heads and an immediate variation of fixed-head flow (spring discharge) in order to balance the water deficit. The resulting pipe flow influences the head gradient of the pipes and, consequently, influences water transfer from the matrix (see Eq. 5) that, in turn, affects matrix heads. The efficiency of this process increases with decreasing hydraulic capacity of pipes (smaller diameter and/or larger roughness) and increasing pipe flows. For the investigated setup, only the conduit diameter of 0.5 m was found to be hydraulically limiting resulting in noticeable head loss along the conduit and, therefore, clearly marked drawdown at the pumping well (Fig. 6). However, the increased water transfer between matrix and conduit induced by conduit drawdown is much smaller than the inflow from the fixed head. Water abstraction from the pipe with 2.5 m diameter does not result in notable conduit drawdown because conduit hydraulics are not limiting. Consequently, water transfer between matrix and pipes is not affected by water abstraction and basically unhampered water inflow through the fixed-head boundary occurs (Fig. 6).
After water abstraction is stopped in period 3, conduit heads immediately rise up to pre-pumping values (Fig. 6). If matrix heads were decreased by preceding water abstraction, i.e., increased water transfer due to conduit head drawdown, actual water transfer will react and matrix heads will return to initial values. This process is reflected by a characteristic delayed rerise of matrix transfer (see Fig. 6 left for the 0.5 m diameter conduit and an abstraction rate of 1.0 m 3 s −1 ).
In summary, this model setup highlights the necessity of a constraining boundary condition for the karst spring to simulate karst water abstraction from conduits. Further, the immediate reaction of conduit drawdown to the onset of water abstraction indicates the necessity of fast storage consideration.

Model enhancement with CADS and the FHLQ boundary condition
The basic model setup (previous section) results in immediate drawdown due to water abstraction because steady-state pipe flow equations do not consider conduit storage. For that reason, the initial model is enhanced by adding CAD storage. The CADS width W CADS is set to 0.25, 0.50, and 1.00 m, respectively, for three different model realizations.
An additional model run with W CADS = 0.00 m is performed for comparison. As previously discussed, the basic model setup demonstrated that in cases of unlimited conduit hydraulics water inflow through the fixed-head boundary dominates (Fig. 6). The fixed-head limited-flow (FHLQ) boundary condition can constrain inflow through the fixed-head boundary resulting in limited inflow to the conduit system. Subsequently, water abstraction of 0.30 m 3 s −1 (node 5) is considered and the fixed-head boundary condition of the basic model setup was extended by a flow constraint (LQ) of 0.025 m 3 s −1 water inflow (node 6, compare Fig. 3). Consequently, the water deficit of 0.10 m 3 s −1 (0.20 m 3 s −1 direct and diffuse recharge minus 0.30 m 3 s −1 water abstraction) is balanced by spring inflow not exceeding 0.025 m 3 s −1 (25 % of the deficit) and additional flow from the matrix continuum via water transfer and from the CAD storage of, in total, 0.075 m 3 s −1 (75 % of the deficit  Bourdet, 1989) s = ∂s ∂ ln t .
Results are presented in Fig. 7. A water deficit occurs with onset of water abstraction (period 2, day 2). Subsequently, the initially fixed-head boundary condition (node 6) switches to limited flow where inflow is restricted to 0.025 m 3 s −1 (Fig. 7a). Without CADS, the model balances the water deficit by instantaneously increased matrix transfer that originates by an instantaneous drop of conduit heads (Fig. 7a, b, see also Eq. 5). If CAD storage is considered, the water deficit is balanced by matrix transfer and CADS flow whereas CADS flow decreases with ongoing time, matrix transfer increases with ongoing time and CADS flow dominates early times (Fig. 7b). Accordingly, conduit heads change less abruptly than without CADS. This effect is increased with increasing CADS width, i.e., more CADS results in less matrix transfer and more damped conduit head drawdown along time. Finally, after a quasi-steady state is reached, matrix transfer is similar to model runs without CADS (Fig. 7a, b). The log-log diagnostic plots (Fig. 7c) clearly indicate the impact of CADS, which acts as karst conduit storage (similar to well bore storage). The presence of this fast storage, indicated by the unit slope of conduit drawdown and drawdown derivative, is increased with increasing CADS (parameterized through W CADS , Fig. 7c). The model without CADS does not show any karst conduit storage in the loglog diagnostic plot.
Water transfer between matrix and conduits affects matrix heads. Consequently, matrix head drawdown is induced by karst water abstraction from conduits. The more water transfer the more matrix drawdown occurs (Fig. 7d), whereas the absolute matrix drawdown is also dependent on matrix parameters (conductivity and storage). It may be noted that initial matrix heads do not depend on W CADS as this parameter does not affect initial conduit heads but only alters the amount of water stored in the CADS. During pumping, however, drawdown in the matrix is slowed down for increased W CADS because more water is abstracted from the CADS in this case.
Subsequently, after water abstraction is stopped (period 3, day 5), the model without CADS does not consider any water deficit. Consequently, the limited-flow boundary at the outlet node 6 switches back to a fixed head and conduit heads recover immediately to pre-pumping values. Because matrix heads are still depressed (Fig. 7d), water transfer from the matrix to the conduits needs some time to recover to the initial value (Fig. 7b, Eq. 5), and spring flow (negative fixedhead flow at node 6) accordingly returns to initial values (Fig. 7a, b). If CADS is considered, the water deficit is still existent because CADS storage needs to be refilled (Fig. 7b, negative CADS flow indicate refilling). Therefore, the constrained boundary condition (node 6) is still active and conduit drawdown recovers with ongoing time in parallel with refilling the CADS (Fig. 7a, b). If CADS is refilled, the water deficit is no longer existent and the limited flow boundary at node 6 switches to fixed head, while matrix heads, and accordingly matrix transfer and spring flow, recover to the initial values. Again, these effects are increased with increasing CADS (Fig. 7a, b).
So far, the analysis demonstrates that both matrix transfer and CADS act in parallel. The sensitivity of model results on CADS is previously investigated. Matrix transfer can be controlled by the transfer coefficient α ex (Eq. 5). Consequently, the sensitivity of model results on this parameter is subsequently investigated. For that reason, an initial model with W CADS = 0.25 m and α ex /l p = 1 × 10 −5 ms −1 is varied by setting α ex /l p to 2×10 −5 ms −1 and 5×10 −6 ms −1 , respectively.
Results are presented in Fig. 8. As already discussed, the onset of water abstraction results in a water deficit that is balanced by matrix transfer and CADS flow (Fig. 8a, b). An increased transfer coefficient results in increased matrix transfer with simultaneously decreased conduit drawdown (Eq. 5) and, therefore, decreased CADS flow (Eqs. 6-8). A comparable behavior occurs for decreased transfer coefficients (decreased matrix transfer, increased conduit drawdown, and increased CADS flow; Fig. 8a, b). The log-log diagnostic plots (Fig. 8c) are reflecting these characteristics with an increasing amount of fast storage, indicated by the unit slope of drawdown and drawdown derivative, for decreasing water transfer coefficients. Further, the variation of the water transfer coefficient is sensitive to initial matrix heads and matrix drawdowns (Fig. 8d). In particular, initial matrix heads are increased for smaller values of α ex /l p , which represent higher hydraulic resistances to matrix-conduit water transfer and, therefore, correspond to larger differences between matrix and conduit heads. In addition, it may be noted that initial matrix heads are again independent of changes in W CADS (cf. Fig. 7d). Matrix drawdown, however, is affected by both parameters. First, increasing W CADS from 0 to 0.25 m (with α ex /l p = 1 × 10 −5 ms −1 in both cases) reduces the drawdown as explained above. Second, variations of α ex /l p alone illustrate that matrix drawdown is lowered for smaller values of this parameter which are also responsible for higher hydraulic resistances to matrix-conduit water transfer, i.e., an increased amount of water is abstracted from CADS during pumping.
This simplified model study demonstrates the importance of the CADS concept: without CADS, any water deficit (for abstraction scenarios) that is not covered by recharge or external boundary condition, e.g., a fixed head, will result in an immediate change of water transfer between conduits and matrix continuum that is directly associated with an immediate variation of conduit heads. CADS provides water immediately (for a water deficit) and, therefore, can dampen the conduit and matrix head variation. It can be concluded that CFPM1 with CADS is able to reproduce the characteristic damped-drawdown behavior within conduits in cases of short-and long-term water abstraction (compare Fig. 1). Log-log diagnostic plots for conduit drawdown and drawdown derivative further approve the existence of fast conduit storage. Overall, CFPM1 with CADS creates testable results.

Cent Fonts case study
A highly idealized representation of the Cent Fonts field situation described by Maréchal et al. (2008) is created to provide an application outlook for using CFPM1 with CADS and FHLQ to represent karst water abstraction scenarios. In doing so, the case study is meant to demonstrate the ability of CFPM1 to reproduce field observations. Data and parameters used for this idealized model are, in most instances, from Maréchal et al. (2008). Several scenarios are performed to investigate parameter sensitivities in response to the onset of water abstraction.

Model setup
The basic features of the study area ( Fig. 1) are conceptualized for modeling purposes as shown in Fig. 9. Cauchy boundary conditions are applied at the north and south borders of the model grid to represent rivers in the catchment area (Fig. 9). The head-dependent water transfer between matrix and rivers is approximated by the MODFLOW River Package with a riverbed conductance of 100 m 2 s −1 . All other lateral outer boundaries of the matrix continuum are of Neumann type (no flow). A uniform diffuse areal recharge of 6.34 × 10 −9 ms −1 (200 mm a −1 ) is applied. The matrix hydraulic conductivity K m is set to 9.00 × 10 −6 ms −1 (to obtain adequate matrix hydraulic heads) and matrix storage is S m = 0.007 (Maréchal et al., 2008). The matrix continuum is discretized by 85 rows and 35 columns with cell lengths and widths equal to 100 m x 100 m. Vertically, the model domain is represented by one unconfined layer with top = 250 m a.s.l. and bottom = −150 m a.s.l.
Highly conductive karst features are represented by one central conduit from north to south, which is subdivided into 90 tubes (each approximately 100 m long) and 91 nodes. CADS is implemented with a width of W CADS = 0.21 m resulting in a storage area W CADS × l CADS of ∼ 1900 m 2 (Maréchal et al., 2008). Conduit node elevation is assumed Hydrol. Earth Syst. Sci., 18, 227-241  at 0 m a.s.l. The conduit diameter is estimated from spring response analysis, according to the concept of Ashton (1966), to be 3.5 m (Birk and Geyer, 2006). Pipe roughness height is set to 0.01 m. Water transfer between matrix and conduit is realized by setting α ex for each node to 4.5 × 10 −5 ms −1 . The water transfer coefficient α ex is doubled in node 1 and node 91 to represent the coupling between river and conduit. The karst spring in the south (node 91) is implemented by an FHLQ boundary condition with fixed head at 76.9 m a.s.l. and inflow limited to 0.03 m 3 s −1 (Maréchal et al., 2008). Water abstraction is realized by pumping from node 87 at 0.4 m 3 s −1 (Fig. 9). Two different time periods are considered: (1) initial period (steady-state) until day 6 and (2) pumping from day 6 to day 38. Beyond the basic model, CADS and conduit-matrix coupling are varied to obtain first insights into sensitivities. Therefore, the CADS width W CADS is set to 0.05 and 0.50 m (basic model 0.21 m), and the water transfer coefficient α ex is varied as 4.0 × 10 −5 ms −1 and 5.0 × 10 −5 ms −1 (basic model 4.5 × 10 −5 ms −1 ).

Fig. 10.
Simulation results for the large-scale pump test scenario: (a) flow terms for variable CADS; (b) flow terms for variable matrixconduit transfer; (c) log-log diagnostic plot for conduit drawdown and drawdown derivatives at node 5; (d) initial matrix head and matrix drawdown at day 38 along the cross section A-A' (Fig. 9), note that initial heads for models with varying W CADS are the same.

Results for the idealized model -initial run and parameter sensitivity
Figure 10 presents flow terms along time, log-log diagnostic plots, and the behavior of matrix heads. In general, water abstraction from the conduit produces a relatively constant drawdown in the pumping well. During the beginning of drawdown formation, CADS flow significantly contributes to balance the pumping induced water deficit. This results in smoother conduit drawdown without an immediate head drop. With ongoing time, CADS flow decreases and matrix transfer increases (Fig. 10a, b). Consequently, CADS influences conduit drawdown and drawdown derivatives for early times much more than matrix transfer. The fast storage, provided by CADS, is reflected by the unit slopes of drawdown and drawdown derivative in the log-log diagnostic plots, too (Fig. 10c). The head gradient within the matrix along cross section A-A' (see Fig. 9) is moderate and the matrix drawdown is less than in the conduit (Fig. 10d). The general behavior of parameter variation with respect to CADS (W CADS ) and matrix transfer (α ex ) is as previously discussed for the simple test model (Sect. 3.2). However, in the context of this application outlook with parameters according to a real situation, it is obvious that the short-term system reaction on hydraulic stress is much more sensitive to CADS (Fig. 10a, b). The smaller the CADS the faster the conduit reaction and the faster a quasi-steady conduit head is reached (Fig. 10a). Further, the quasi-steady conduit head depends on the conduit-matrix coupling, here varied via the transfer coefficient α ex (Fig. 10b). In fact, drawdowns in the conduit pumping well are reduced with better coupling (increased α ex ) because the necessary head difference between matrix and conduit to result in a certain water transfer is reduced (see also Eq. 5). On the contrary, smaller water transfer coefficients result in enhanced conduit drawdown (Fig. 10b).
The initial matrix head distribution is sensitive to the transfer coefficient because this parameter regulates flow and head difference between matrix and conduits (Eq. 5; Fig. 10d). Further, the current system understanding indicates that the initial matrix head distribution depends on the spatial distribution of the conduit network and the transfer coefficients. Due to our conceptual model, the matrix is mainly drained by conduits and, therefore, variations of α ex are strongly affecting matrix heads. Initial matrix heads are not influenced by CADS because steady-state situations do not account for storage. However, under the used parameters, matrix drawdown is sensitive to CADS and clearly decreases with increasing CADS (Fig. 10d) because matrix-conduit transfer and conduit drawdown is decreased (Fig. 10a). Matrix drawdown is comparatively less sensitive to matrix-conduit transfer (Fig. 10d).
Hydrol. Earth Syst. Sci., 18, 227-241, 2014 www.hydrol-earth-syst-sci.net/18/227/2014/ Water budget terms are given in Table 1. Hence, for the initial model about 10 % of water pumped during period 2 comes from CADS and about 38 % comes from matrix storage. Further budget terms for the parameter study underline the already discussed model behavior. Accordingly, the amount of pumped water coming from the matrix storage decreases as CADS increases (Table 1), as expected. Furthermore, the conduit-matrix coupling (α ex ) does not significantly affect the distribution of water coming from CAD storage and from matrix storage, see Table 1.
In principle, CFPM1 with the CADS and FHLQ functionality is able to qualitatively reproduce the field situation described by Maréchal et al. (2008) (Fig. 1). It can be concluded that CADS has a strong influence on model reaction to hydraulic stresses like the onset of water abstraction. On the contrary, matrix-conduit transfer is very sensitive to the initial matrix heads. Matrix heads and matrix drawdown vary significantly with distance from the conduit. Consequently, matrix heads seem to be very valuable to estimate the spatial distribution of the conduits.

Comparison with measured data
Subsequently, the model is further adapted according to the situation described by Maréchal et al. (2008) in order to compare model results with field measurements. The previous analysis demonstrates that the system strongly reacts to hydraulic stress. Because the pumping well is directly placed in the highly conductive conduit, pumping rate variations are expected to strongly affect hydraulics. Consequently, measured pumping rates, slightly variable with time (see Fig. 1 right), are considered by CFP as time-dependent input data with a resolution of t = 3600 s. Further, over period 2 (pumping) the diffuse areal groundwater recharge is reduced to 10 % of the initial value (6.34 × 10 −9 ms −1 ) because the field experiment was conducted during a dry period without recharge (Maréchal et al., 2008). Here, the remaining recharge of 6.34 × 10 −10 ms −1 is assumed as background value due to slow draining of the less conductive rock matrix.
Two model setups are automatically calibrated using PEST (Doherty, 2005) whereas K m (matrix hydraulic conductivity, upper and lower boundary 1.00 × 10 −5 ms −1 and 1.00 × 10 −8 ms −1 ), S m (matrix storage, upper and lower boundary 2.00 × 10 −1 and 1.00 × 10 −4 ), and α ex (transfer coefficient, upper and lower boundary 1.00 × 10 −4 m 2 s −1 and 1.00 × 10 −8 m 2 s −1 ) are considered as free parameters. CADS is not considered for calibration in order to reduce the number of free parameters. Rather, CADS is parameterized according to Maréchal et al. (2008). Hence, setup (1) uses W CADS = 0.21 m. Setup (2) is intended to investigate how CFP without CADS can reproduce field observations. Consequently, CADS is deactivated by setting W CADS = 0.00. Calibration considered measured conduit drawdown at the pumping well (node 87) plus matrix drawdown. Because the position of matrix drawdown relative to the conduit is unknown, only a rough estimation of h m = 5 m (over the pumping period) with h m,ini = 110.0 m a.s.l. (Maréchal et al., 2008) is assumed at position M1 with a distance of 1000 m between conduit and matrix observation well (Fig. 9).
Computed flow terms and conduit heads are presented in Fig. 11. In general, both models show similar behavior of computed and measured conduit drawdown. However, the model without CADS fails to reproduce periods with highly variable pumping rate (onset of water abstraction at day 6, interrupted pumping around day 14) due to the missing fast storage. On the contrary, CFP with CADS does represent much better the early stage of the pumping test as well as the interrupt around day 14. Further, it is obvious that pumping rate variations considerably influence flow terms. For the CFP model without CADS, matrix-conduit transfer responds to pumping rate variations. As previously discussed (Sect. 3), this is closely connected with an immediate conduit head variation. If CAD storage is accounted for, mainly CADS flow responds to pumping rate variations whereas matrixconduit transfer is basically unaffected. Consequently, the conduit head drawdown curve is much smoother due to the CADS caused damping.
The resulting log-log diagnostic plots (Fig. 11b) underline the significance of CADS in order to represent fast storage that results in the unit slopes for conduit drawdown and drawdown derivative. Keeping in mind that the aim of the simplified model is to evaluate different model concepts rather than to find a good fit (which strongly depends on assumptions on conduit geometry), the deviation between measured and modeled heads is acceptable for setup (1). Figure 11b also shows that matrix head drawdown computed at observation well M1 (for location see Fig. 9) is qualitatively similar to matrix head drawdowns reported by Maréchal et al. (2008) as indicated in Fig. 1. Therewith, CFP with CADS is able to qualitatively describe the drawdown behavior with time for both conduit and matrix heads. On the contrary, setup (2) is not able to represent the initial phase of the pumping test as well as the reaction of conduit heads on strong variations of the pumping rate because the model lacks CADS (Fig. 11a,  b). This produces the already described instantaneous drop of conduit heads missing the initial period of unit conduit drawdown and drawdown derivative (Fig. 11b).
The matrix heads prior to pumping as well as the drawdown after 38 days of pumping along the A-A' cross section are shown in Fig. 11c. Due to the increased transfer coefficient in setup (2), the initial hydraulic gradient within the matrix is more pronounced. Because setup (2) lacks CADS, the matrix drawdown is increased, too. The following parameters are obtained after calibration: for setup (1) with W CADS = 0.21 m: α ex = 8.72 × 10 −5 m 2 s −1 , K m = 3.00 × 10 −6 ms −1 , S m = 0.0011; for setup (2) with W CADS = 0.00 m: α ex = 2.85 × 10 −4 m 2 s −1 , K m = 1.00 × 10 −6 m s −1 , S m = 0.0007. As the model without CADS needs to infer matrix-conduit transfer to respond to pumping induced hydraulic stress (Fig. 11a), α ex is increased to achieve a better conduit-matrix coupling. Further, K m as well as S m are decreased to increase Hydrol. Earth Syst. Sci., 18, 227-241, 2014 www.hydrol-earth-syst-sci.net/18/227/2014/ the matrix responding behavior of the model without CADS. Therewith, the gradient of matrix heads is steeper and matrix head drawdown is increased (Fig. 11c). Finally, it can be summarized that CFP with CADS is a suitable tool to represent karst catchments under water abstraction (and other hydraulic stresses). Thereby, CADS is necessary to account for the fast storage component associated with conduits that results in non-abrupt drawdown of conduit heads. Further, the case study highlights that matrix heads are sensitive to (a) the transfer coefficient, (b) the distance from the conduit as well as (c) CADS (Fig. 11). Several modifications appear to be possible to achieve further model improvement, for example a more realistic consideration of the model domain geometry, the local geology, and the distribution of karst conduits within the matrix. This may help to further reduce deviations between modeled and measured conduit heads (Fig. 11a). Additional analysis can be applied for further model evaluation like flow dimension analysis as indicated by, for example, Walker et al. (2003), Maréchal et al. (2008), and Cello et al. (2009).

Conclusions
Implementation of conduit-associated drainable storage (CADS) to the existing Conduit Flow Process Mode 1 (CFPM1) combines the conceptual approaches for water storage in karst systems presented by Mangin (1975Mangin ( , 1994 and Drogue (1974Drogue ( , 1992 resulting in a triple porosity system representation (Worthington et al., 2000). Thereby, matrix and fracture porosity are usually merged within one continuum because an REV can be defined. Hydraulic parameters of the REV (fissured/fractured matrix blocks) can be obtained from traditional hydraulic borehole tests (e.g., Geyer et al., 2013). Fast reacting conduit-associated storage is represented by CADS. The newly developed functionality is fully integrated in the CFPM1 flow subroutines and requires only the storage width as an additional model parameter. The CADS volume (parameterized by the CADS width W CADS ) is a calibration parameter with a physical background and can be obtained via transient model calibration, for example, from the reaction of conduit heads on hydraulic stresses like start of pumping, stop of pumping, or strong recharge signals directly routed in the conduits. As CADS volume is conceptually independent from the conduit volume, the chosen approach allows a rather flexible treatment of water storage or release behavior linked to the highly permeable flow system.
The CADS approach was evaluated in several pumping test scenarios to investigate the effect of storage properties and boundary conditions on karst hydraulics. Simulation results show that associated conduit storage plays a major role during the early time periods of water abstraction from karst systems, even though the majority of total aquifer storage is provided by the matrix storage. This is primarily due to the fact that only a small amount of water from matrix stor-age is provided during the early stages of water abstraction through conduit-matrix coupling. The newly implemented CADS flattens the drawdown curve from the beginning of water abstraction from the conduit system because of immediate water inflow from CAD storage. Depending on the model setup, strong water abstraction within highly permeable structures, i.e., conduits, can result in unhampered water inflow through fixed-head boundaries connected to the pipe network. This effect leads to minor drawdown during water abstraction even at high pumping rates and consequently an insignificant contribution from the CADS. This condition can be relevant, for example, for streams that are hydraulically perfectly connected to karst conduit systems. Therefore, the simulation of pumping tests with CFPM1 can require the additional implementation of a fixed-head limited-flow (FHLQ) boundary condition, which constrains inflow for constant head boundaries. For these scenarios, water deficit resulting from water abstraction from the conduit system is balanced out by contributions from the CADS and the matrix storage.
CADS is assumed to be uniform (width is constant with depth). Furthermore, groundwater flow within CADS (e.g., horizontal or vertical flow driven by hydraulic gradients within the pipe) is not represented. Wellbore storage also is not considered by CADS, however, the parameters of CADS may be adjusted during calibration to account for wellbore storage. CADS and CFPM1 updates that overcome these limitations will eventually be available under the subsequently provided internet domain.
CADS and the FHLQ boundary were further evaluated to simulate a large-scale field pumping test reported by Maréchal et al. (2008). Dimensions and hydraulic model parameters were set in the range of observed field values. Even though the geometry of the karst aquifer was highly idealized, the model was able to qualitatively reproduce the overall drawdown curves observed in the pumping well and the observation wells. Initial comparison of computed and measured drawdown demonstrates the necessity of fast storage associated with conduits. Ongoing work will evaluate the large-scale pumping test with CFPM1 and CADS by using an adequate hydrogeological representation of the catchment and further diagnostic tools like drawdown derivatives and flow dimension analysis. Further work will focus on systematic-type curve analyses to evaluate pumping test responses under different complex modeling setups. The newly developed CADS package and FHLQ boundary can be used for evaluation of large-scale karst aquifer hydraulic tests with relatively modest input data requirements. Further, even more ambitious studies concerning, for instance, spatially or temporally variable recharge or hydraulic and hydrochemical responses of karst springs as outlined by Hartmann et al. (2012Hartmann et al. ( , 2013 or Long and Mahler (2013) appear to be feasible in the long run.