EXAMINING A COUPLED CONTINUUM PIPE-FLOW MODEL FOR GROUNDWATER FLOW AND SOLUTE TRANSPORT IN A KARST AQUIFER PREIZKUS MODELA SKLOPLJENEGA TOKA V ZVEZNEM POROZNEM SREDSTVU IN V CEVEH ZA TOK PODZEMNE VODE TER PRENOSA SNOVI V KRAŠKEM VODONOSNIKU Bill x. HU1 Abstract UDC 556.33 Bill X. Hu: Examining a Coupled Continuum Pipe-Flow Model for Groundwater Flow and Solute Transport in a Karst aquifer A coupled continuum pipe-flow (CCPF) model has been developed for groundwater flow and solute transport in a karst aquifer with conduits. Groundwater flow in conduits is simulated through a pipe flow model and flow in fissured matrix rock is described by Darcy's law. Water mass exchange between the two domains is modeled by a first-order exchange rate method. In this study, we investigate mathematical well-posed-ness (mathematical term, which means solution existence and uniqueness) of the CCPF model, develop a finite elementary method to numerically approximate the mathematical model and study the convergence of the numerical method. The study results prove the modeling approach is mathematically well posed and numerically converged. To study the accuracy of the CCPF model, a recently developed Stokes-Darcy (SD) model and CCPF model are compared with laboratory experimental results. It was found that the SD model simulations match well with experimental results, but the CCPF model overestimates the hydraulic head in the matrix, especially around the matrix and conduit interface. The model underestimates solute transport in the conduit and does not capture the plume distribution in the matrix. In comparison with the SD model, the CCPF model requires an additional parameter, the first-order mass exchange rate, and the parameter is normally obtained through inverse method curve fitting. The SD method may provide an approach to directly estimate the parameter value. Keywords: karst aquifer, continuum pipe-flow model, Stokes-Darcy model, mathematical well-posedness, mass exchange rate. Povzetek UDK 556.33 Bill X. Hu: Preizkus modela sklopljenega toka v zveznem poroznem sredstvu in v ceveh za tok podzemne vode ter prenosa snovi v kraškem vodonosniku Za tok podzemne vode in prenos snovi v kraškem vodonosniku s kanali je bil razvit model sklopljenega toka v zveznem poroznem sredstvu in v ceveh (CCPF). Podzemni tok v kanalih je bil simuliran z modelom toka v cevi, tok skozi razpoklinsko matrico pa opisuje Darcyjev zakon. Masna izmenjava vode med dvema domenama je bila modelirana z metodo hitrosti masne izmenjave prvega reda. V tej študiji smo preizkušali matematično dobro postavitev modela CCPF (matematični izraz, ki izraža obstoj rešitve in njeno edinstvenost), razvili dokončno osnovno metodo za numeričnen približek matematičnemu modelu ter proučili ujemanje z numerično metodo. Rezultati študije dokazujejo, da je modelarski pristop dobro postavljen in se numerično ujema. Da bi ugotovili natančnost modela CCPF, smo nedavno razviti Stokes-Darcyjev model (SD) in model CCPF primerjali z rezultati laboratorijskega eksperimenta. Ugotovljeno je bilo, da se simulacije modela SD dobro ujemajo z rezultati poskusa in da model CCPF precenjuje piezometrično višino v matrici, še posebej na območju meje med kanalom in matrico. Slednji model tudi podcenjuje prenos snovi v kanalu in se ne ujema dobro s porazdelitvijo koncentracij raztopine v matrici. V primerjavi z modelom SD zahteva CCPF dodatni parameter hitrost masne izmenjave prvega reda. Ta parameter pa navadno pridobimo z inverzno metodo prileganja krivulje. Metoda SD lahko omogoči pristop k neposrednemu izračunu vrednosti tega parametra. Ključne besede: kraški vodonosnik, sklopljen tok v zveznem poroznem sredstvu in v ceveh, Stokes-Darcyjev model, matematično dobra postavitev modela, hitrost masne izmenjave. 1 Department of Earth, Ocean and Atmospheric Sciences, Florida State University, Tallahassee, FL 32306, USA, e-maü: hu@gly.fsu.edu Received/Prejeto: 16.10.2009 INTRODUCTION Karst aquifers supply about 25% of the world population with water (Ford & Williams 1989), including nearly all the drinking water to certain regions, e.g., 90% of Florida's drinking water (Paulson et al. 1990). The presence of conduits is unique to karst aquifers. Phreatic conduits are connected, water-saturated tubes located below the water table. These conduits largely control groundwater flow and contaminant transport within the aquifer (Katz et al. 1998). The water-saturated matrix rock surrounding the conduits comprises most of the aquifer volume and hence most of the storage (Worthington 2007). It has been observed in field sites that during a high-flow event, the water pressure in the conduits is larger than that in the ambient matrix so that conduit-borne contaminants can be driven into the matrix. During a low-flow event, the pressure differential reverses and contaminants sequestered in the matrix can be released into the free flow in the conduits and exit through springs into surface water systems (e.g., Martin & Dean 1999, 2001; Li et al. 2008). This retention and release phenomenon induces an environmental issue in that sequestered contaminants may influence the quality of underground water sources for a long time and thus significantly decrease quality water availability. The dual character of a karst flow system is widely recognized and stems from the existence of different porosities within a karst aquifer (Ford 2003; Worthington 2003). The porosity difference determines the type of flow prevailing in the aquifer (Ford & William 1989; Bauer et al. 2000). Similar to the dual-porosity/permeability model widely used for fractured media (e.g., Gerke & van Genuchten 1993a, b), the coupled continuum pipe-flow (CCPF) model has been proposed to describe the flow and solute transport in karst aquifers (Chen & Bian 1988; MacQuarrie & Sidicky 1996; Kiraly 1998; Bauer et al. 2000, 2003; Taylor & Greene 2001; Birk et al. 2003). The CCPF model is a dual flow system consisting of a matrix representing the bulk mass of permeable limestone and a conduit system representing the karst conduit network. Flow exchange between the two systems is controlled by differences in hydraulic heads as well as the hydraulic conductivity and the geometric setting. In the CCPF model, groundwater flow in the matrix is described by Darcy's law, and flow in the conduit is modeled by a pipe-flow model. The water mass exchange flow rate between the two systems, qex, is described by a firstorder mass exchange model. The exchange flow rate is assumed to be linearly proportional to the head difference between the two systems (Barenblatt et al. 1960; Cao et al. 1988; Sauter 1992; Teutsch 1989). The exchange rate coefficient is a lumped parameter. Its value will depend on many factors including: hydraulic conductivity in the matrix, the exchange surface between the conduit and matrix, and conduit geometry (Barenblatt et al. 1960; Liedl et al. 2003). The value of the exchange rate parameter is not usually obtained from measurements but rather through curve-fitting. Based on the CCPF model, a new numerical method has been developed and has become part of the new MODFLOW software (Shoemaker et al. 2008). However, the suitability and validity of the CCPF model as a model for groundwater flow in a karst aquifer, especially for the flow exchange between matrix and conduits, has not been well studied. In addition, determination of the value of the exchange rate parameter is also an issue that needs attention. Recently, Faulkner et al. (2009) developed a numerical model to simulate groundwater flow and solute transport in a karst aquifer based on a dual-regional conceptual model. The karst aquifer is divided into two regions, the limestone matrix and the conduits. Ground-water flow in the matrix is still assumed to satisfy Darcy's law, but conduit flow is described by the Stokes equations. The renowned empirical Beavers-Joseph interface condition (Beavers & Joseph 1967) and its simplified version, the Beavers-Joseph-Saffman condition (Saffman 1971), are used to model the flow exchange at the interface between the matrix and the conduit. Faulkner et al. (2009) also conducted a sandbox experiment to simulate groundwater flow and solute transport in a karst aquifer with a single conduit. The experiment results match well with the Stokes-Darcy simulations. In this study, we will first introduce the CCPF model, mathematically investigate its well-posedness (mathematical term, means solution existence and uniqueness) and regularity, develop a numerical finite elementary method for a special case, and evaluate the numerical convergence of the model. The CCPF and SD numerical simulations will then be compared with the experimental results by Faulkner et al. (2009). CONTINUUM PIPE FLOW MODEL Conventionally, hydraulic interaction between the fissured system and the conduit network has been simulated by employing three types of modeling approaches (Liedl et al. 2003; Birk et al. 2003; Faulkner et al. 2009). First, multiple sets of fractures may be coupled to each other in order to represent the different hydraulic properties of the fissured and the conduit system. Alternatively, double-continuum or multi-continuum models may be applied where the cross-flow between the flow systems depends on the corresponding pressure differences via linear exchange terms. As a third approach, discrete networks of flow paths may be coupled to a continuum in order to model the dualistic flow patterns in karst systems. The last one of the three approaches is the most popular among hydrologic studies, called the coupled continuum pipe flow model. As the name implies, this is a coupled system consisting of a continuum, the matrix, with a pipe flow conduit imbedded inside. It has been used successfully to study groundwater flow in karst aquifers and the genesis of karst aquifers (e.g., Bauer et al. 2003; Liedl et al. 2003; Birk et al. 2003). In this section, we will develop a three-dimensional version of this model. For conduit flow, the discharge can be related to the head difference in the tube by applying the Darcy-Weis-bach equation (Bobok 1993): 64v dhc _ ^ "" dr ~ 2dg (1) where t is the tangential unit vector along the 1D pipe conduit, d is the pipe diameter, u=4Q/{7id^) is the average velocity, and g is the earth's gravitational acceleration. Here, Q is the total discharge in the pipe. The friction factor, X, depends on the velocity in the pipe via the Reynolds number Rg= ud/v, with v the kinematic viscosity of water. For a low flow velocity, laminar flow is assumed and the Hagen Poiseuille equation can be applied. The friction factor for laminar flow is calculated as A =-. Plugging the above relations into du (1), we have Q = ^ 128v dr Since dQ/dr = -/„,, we have a nd'^g dhc dr^ 128v dr^ (2) The fissured matrix and the conduit system are coupled at common nodes by a quasi-steady-state exchange term (Barenblatt et al. 1960), fm=<^exiK -hJ (3) Where aex > 0 is the water mass exchange coefficient. This equation indicates that the process in the conduit is controlled by groundwater flow in the fissured matrix. The exchange coefficient is the key parameter, into which, unfortunately, most of the uncertainties are lumped. A detailed discussion of the conventional wisdom in determining aex will be presented later. In the steady-state case, we have the following linear system, which is essentially two coupled Poisson equations in different domains. inn. (4) where £)_"** ž , jand g denote the external forces, and 128 V 0„; denote the regions for the fissured continuum and pipe conduit, respectively. WELL-POSEDNESS AND REGULARITY The time dependent version of equation (4) has been used to study karst aquifer genesis (Bauer et al. 2000, 2003; Bauer et al. 2003). The model was originally used for cou- pling non-linear Richard's equation and pipe flow, which studies variably saturated media (MacQuarrie & Sudicky 1996). The mathematical model is solved by the Carbon- ate Aquifer Void Evolution (CAVE) program (Clemens et al. 1996). CAVE solves the flow in the fissured matrix by a finite difference scheme using MODFLOW (Har-baugh 2005) and the conduit flow by a non-linear finite difference discretization (The non-linearity only occurs for the turbulent case. Since we are only considering laminar flow in the conduit, a non-linear solver is not necessary). The coupling exchange condition is reached by iteratively solving equations in the matrix and conduit regions until convergence. The exchange of fluid is only allowed at discrete nodal points. No rigorous mathematical theory was presented therein to show well-posedness or guarantee convergence. Here, the problem is simplified. We study a stationary case with linear conduit flow. It is then worthwhile to systematically present the mathematical theory behind the model. We assume the following geometrical setting for our problem. The fissure continuum is assumed to occupy the unit square = (0,1) x (-1/2, 1/2) C R2 and the one dimensional conduit pipe lies in the middle Ü.C = (0,1) x {0}, i.e., a horizontal straight pipe in the middle of the unit square matrix. Then equation (4) becomes -V{KWhJ=a,^{h,-hJS{y) + f in Cl^ = in 5 ax ax where 8(y) is the Dirac delta function in y. In addition, we use the following Dirichlet boundary conditions for both domains. K=gm on an^ K =Sc on ^^c where gm and g, are given functions. We remark that we may also consider a three dimensional fissured matrix coupled with two dimensional free flow plane or one dimensional free flow pipe. The well-posedness proof would remain the same, while the regularity and numerical analysis could be treated in a very similar fashion, although the results slightly differ. In the following analysis, we only consider the homogeneous case, i.e., we assume that gm = 0 and g, = 0. The non-homogeneous boundary case can be converted to the homogeneous case through a standard procedure. First, we formulate the variational form for equation (6). To this end, we define a bilinear from a(- ,•) on H:=Ho(n,„)xHj(n,.) as follows. For h = (hm, h,) and v = (v^, v,) in H, a(h,v):=f^KVh^ix, y)-Vv^ix,y)dxdy + +jIdW, {x)V, ix)dx 7 -«exfo^K (x>0) - h^(x))v^(x.O)dx + a^f^(hc (x.0) --h^(x})Vc(x)dx We say that h GH is a weak solution of equation (6) if a(h,v) = (f,v„,) + (g,v,) V v(EH (8) If we assume that ^eH-1(Om) and geH-1(0,), then weak solution of (6) exists and is unique. We will prove that the bilinear form is continuous and coercive in H. The existence and uniqueness then follow from Lax-Milgram theorem. The continuity follows directly from the trace theorem. For h (E H, fl(h,h) = KVh^ U y) ■ Vh^ix, y)dxdy + + flD{h[ {x,0)-h„ix))h^{x,0)dx + oc^f^(h,(x,0) -h^(x))h,(x)dx>fn KVh^ ix, y) • ■Vh^ix.y)dxdy + f^DiK {x)fdx which implies the coercivity. Let h = (hm, hc) be the weak solution of (6). Assume that and ^EL^CfiJ. Then where Q==f(h,{x,0)--h^(x)Mx,0)dx -i-E Thus S{y){h,.^ (n^). By equation (6) we conclude that (n,„). FINITE ELEMENT APPROXIMATIONS AND CONVERGENCE Now we construct finite element approximations for the weak solution of equation (6). Let and Oj be a quasiuniform triangulation of 0m and 0c, respectively. We construct the conforming finite element spaces Hmh and Hch as piecewise continuous quadratic function spaces of Hl(0m) and H^(0c), respectively, where h is the mesh size. The finite element approximation for the weak solution h of equation (6) is to seek hh=(hh„, hC)GHm x Hh such that To study the convergence of the numerical solution, we set all parameters to one. We adjust the forcing f and g in matrix and conduit, respectively, such that the following solution is exact. lij = 2sin(27Ec) in Q^ h,„ = sin(27rjc) in (0,1) x (-1 / 2,0) D D ,„ (10) hm={-0Cexy/K + l)sm{2nx) in (0,1) x [0,1/2] DQ^ Tab. 1: Convergence rate for steady state CCPF. (9) h Il hc - hh ||0 || hm - hm ||0 | hc - hh |0 | hm - hm |0 2-3 3.047E-2 1.644E-2 7.888E-1 6.906E-2 2-4 3.907E-3 2.092E-3 2.025E-1 1.517E-2 2-5 4.915E-4 2.628E-4 5.096E-2 3.582E-3 2-6 6.150E-5 3.290E-5 1.276E-2 8.818E-4 2-7 7.694E-6 4.114E-6 3.191E-3 2.198E-4 rate of conv. 2.989 2.992 1.987 2.070 We have the following convergence rate as is summarized in Tab. 1. tte better-than-predicted convergence rate is because we have the edges of the elements in 0mh to lie exactly on 0ch, i.e., we do not have any of 0hc to intersect the interior of the elements in 0mh . Furthermore, it is speculated that h has piece-wise higher regularity. tte velocity determined from the CCPF system is used in the governing equation for the solute evolution in the matrix: By the regularity result and a standard argument for finite element approximations, we have the following error estimate. Equation (9) admits a unique solution hh=(hm, hh)GHm x HJ. Moreover, for 0 < £ < -j, there exists a constant C(e), which is independent of h such that h-h^ dt (11) where Cm denotes the solute concentration in matrix and D is dispersion coefficient. INTRODUCTION OF STOKES AND DARCY MODEL Faulkner et al. (2009) developed a SD model to simulate groundwater flow and solute transport in a karst aquifer with a conduit in matrix. Here we would like to briefly introduce the method for the completeness of the paper. If interested, the reader could find the detailed description of the method from Faulkner et al. (2009). ^e flow in the matrix 0m is governed by the Darcy system hJx,t = Q)=h, in 0 (12) where qm denotes the specific discharge, which can be expressed as qm = nvm, where vm denotes the seepage velocity and n the effective porosity, S the storage coefficient, K the hydraulic conductivity tensor, and fm represents sink/source term. tte hydraulic head hm is defined b^ ^^ ^ + where pm denotes the pressure, r o p the water density, g the gravity acceleration, and z the position head. By substituting the second equation in (12) into the first one, we obtain the equation that governs the change of the hydraulic head: dt + y.(-KWhJ=f,„ in 0„ (13) tte boundary conditions for (13) are the Dirichlet boundary condition hm=H(x) along rg, where h(x) is known from measurements, and the homogeneous Neumann boundary condition (KVhm) • n=0 along r0 that represents a no-flow boundary condition at the artificial boundary of the aquifer. In the conduit domain, Oc, the other domain of the problem, as demonstrated by Faulkner et al. (2009), the Navier-Stokes equations govern the free flow, f-v-r.x V-v, =0 in Oc (14) where vc denotes the fluid velocity, T(v, p)=-pci+2vD(v) the stress tensor, p^ the kinetic fluid pressure, D(v)=^(Vv+(Vv)T) the deformation tensor, v the kinetic viscosity of the fluid, and fc a general body forcing term. If Re is small, the advective term (v •V)v is negligible and equation (14) could be simplified. Faulkner et al. (2009) developed suitable boundary conditions for the system. At the sinkhole and the spring, nonhomogeneous Dirichlet boundary conditions are applied to specify the inflow and outflow velocities, respectively. Specifically, VcXn=0 and Vc-n=ysi(t)nsi(x)=fsi on r^; VcXn=0 and v^ •n=ysp(t)nsp(x)=fsp on r^^ (15a) (15b) where ysi, y,^, nsi and n5p are given functions defined at the spring, rjp, and sinkhole, r,,. ^ese boundary data are obtained from field measurements. tte right column corresponds to the Dirichlet normal velocity. tte left column corresponds to the Dirichlet tangential velocity because nxvcxn=vc-(vc •n)n. In addition to the boundary conditions imposed along the boundaries of the matrix and conduit, the following interface boundary conditions are used by Faulkner et al. (2009) to couple the solutions in the two nonoverlapping yet neighboring domains: -ncmT(Vc.p)ncm =g(h^-z) yj trace (k) onI,,„ (16) Where rcm denotes the conduit-matrix interface, T represents the local tangent plane to rcm, ncm denotes the unit normal vector to rcm pointing from the conduit to the matrix, a denotes a constant and k represents the permeability, which has the following relation with the hydraulic conductivity, K=k g/v. It should be noticed that k and K differ by a factor of a constant scalar for a certain type of fluid. ttus, all assumptions on K such as the symmetric positive definiteness also carry over to k. tte velocity determined from the coupled Stokes-Darcy system is used in the governing equation for the tracer evolution in the matrix, equation (11). COMPARISON BETWEEN MODELING SIMULATIONS WITH LABORATORY EXPERIMENTAL RESULTS In our previous study (Faulkner et al. 2009), we used a laboratory analog experiment to simulate groundwater flow and solute transport in a karst aquifer with one conduit buried adjacent to the matrix. tte experiment mainly focuses on the water and solute exchanges between the matrix and conduit. tte study results will be used as the benchmark for numerical study. tte experimental facility and procedure are described in Faulkner et al. (2009). Here, we just use the results to compare with the numerical simulations. Fig. 1 shows the experimental hydraulic head distribution and simulation results by SD (Faulkner et al. Fig. 1: Experimental (top), SD simulated (middle) and CCPF simulated (bottom) head distributions in the matrix, where experimental and SD results are from faulkner et al. (2009), and CCPF results are new. The head unit is centimeter. fig. 2: Experimental (left), SD simulated (middle) and CCPF simulated (right) results for the solute concentration distribution in the matrix at various time instants; top to bottom: t=32.5s, t=62.5s, t=92.5s and t=122.5s, where experimental and SD results are from faulkner et al. (2009), and CCPF results are new. ^e simulated concentration values change from 1.0, the deep purple, to 0.0, white. 2009) and CCPF models, respectively. The hydraulic heads from the laboratory measurements and simulations by SD and CCPF at measurement points are listed in Tab. 2. From the Fig. 2 and Tab. 2, we can tell that the SD simulations are close to the experimental results; while the CCPF modeling generally overestimates the hydraulic heads in the matrix, especially at the boundary between the matrix and the conduit. Fig. 2 presents the experimental results of the dye distributions at several time slots and simulations at the same time slots by the two models. The SD simulations are very similar to the experimental results, but the CCPF modeling results are quite different. First, the pipe-flow modeling underestimates dye front movement in the conduit as well as dye exchange at the front between the two domains. Second, the modeled plume distribution in the matrix is convex shaped and does not capture characteristics of the plume distribution in the experiment (a broad U-shaped plume with two humps at either end caused by end point interface). Tab. 2: Comparison of experimental hydraulic head results with SD and CCPF simulations at measurement points. location (cm) Lab results SD Model Diff Between SD CCPF Diff Between CCPF rOinl (x, y) (cm) (cm) and Lab (cm) Model (cm) and Lab (cm) 1 (55.3,13.2) 23.09 23.10 0.01 23.10 0.01 2 (47.0,6.5) 25.75 28.46 2.71 27.39 1.64 3 (47.0,13.2) 25.48 26.39 0.91 26.32 0.84 4 (47.0,19.9) 25.21 25.71 0.50 25.85 0.64 5 (37.8,6.5) 30.06 32.12 2.06 32.70 2.64 6 (37.8,13.2) 29.26 30.03 0.77 30.61 1.35 7 (37.8,19.9) 29.61 29.05 -0.56 29.63 0.02 8 (28.6,6.5) 32.29 32.96 0.67 35.14 2.85 9 (28.6,13.2) 30.95 31.58 0.63 32.94 1.99 10 (28.6,19.9) 30.15 30.81 0.66 31.84 1.69 11 (19.4,6.5) 32.67 33.03 0.36 35.15 2.48 12 (19.4,13.2) 31.13 31.85 0.72 33.17 2.04 13 (19.4,19.9) 29.75 31.26 1.51 32.24 2.49 14 (10.2,6.5) 31.17 31.81 0.64 32.44 1.27 15 (10.2,13.2) 30.56 31.02 0.46 31.54 0.98 16 (10.2,19.9) 30.11 30.73 0.62 31.14 1.03 17 (2.5,13.2) 30.20 30.21 0.01 30.21 0.01 Summation of errors T 12.68 23.97 Summation of absolute errors 13.80 23.97 SENSITIVITY STUDY ON aex AND a As is pointed out and repeatedly verified in previous studies (Bauer et al. 2000, 2003; Birk et al. 2003), coupled continuum pipe flow models are sensitive to the choice of the first order exchange parameter aex. tte breakthrough time of conduit genesis may vary among several magnitudes as aex varies in the range of hydraulic conductivity. In this study, we use the discharge and discharge boundary condition for the CCPF. For this pair of boundary condition, the CCPF model performs rela-lively well. Quantities such as the exchange of fluid along the interface and the discharge in the conduit are sensitive to the parameter choice. In the laboratory experiment, the hydraulic conductivity of the glass beads is, K~7.4x10-4 (m/s). So we set aex in the range of [10-1, 10-4] (m/s). It is observable from Fig. 3 that the exchange flow is sensitive to the choice of aex. For a comparison with the SD system, the results are summarized in Fig. 4, where aex ranges between 10-4 to 10-2 (m/s). From the results, we see that, for conventionally suggested small aex, the conduit flow is almost Fig. 3: Variation of normal velocity along the interface as aex in the CCPF simulation. Fig. 4: Comparison between Stokes-Darcy (green) and CCPF (blue) simulation results of exchange flow (left) along interface and conduit discharge (right) under various aex values: a) aex =10-4 m/s; b) aex =10-3 m/s; and c) aex =10-4 m/s. Fig. 5: Sensitivity study on a used in Beavers-Joseph interface condition: variations of normal velocity (left) and conduit discharge (right) with a change from 0.2 to 2.0. uniformly linear and interface fluid exchange is almost a constant, which is far from the SD model results. tte SD model can be regarded as the "correct model" since its validity is justified in Faulkner et al. (2009). ^e high sensitivity at small aex value can be understood in the following way. In the case of ae = 0.0 m/s, Darcy's flow and pipe flow are decoupled and hence there should be a non-zero head difference in general. tte sensitivity equation with respect to aex (which can be derived by formally differentiating the CCPF system with respect to a) takes the same form as CCPF but with the head difference serving as the forcing/source term. Since the hydraulic conductivity and aex are small, the response (the sensitivity) is expected to be large. On the other hand, one may speculate that the Stokes-Darcy system is sensitive to the choice of a in the Beavers-Joseph interface condition as well. However, this is not the case. From Fig. 5, we do not observe a visible difference if we vary the a in the Stokes-Darcy model according to the range suggested by Beavers and Joseph (Beavers & Joseph 1967). Study results indicate that the a in the Stokes-Darcy model with Dirichlet boundary conditions is a "dummy" parameter, and simulation results are very insensitive to its variation. tte parameter value, in practice, can be assumed to be "known". SUMMARY AND CONCLUSIONS ttis study has examined a currently used modeling approach, CCPF, for groundwater flow and solute transport in a karst aquifer with conduits. We have examined whether the modeling system is mathematically well-posedness, and then have developed a finite numerical method to solve the mathematical model and studied the convergence of the numerical method. tte modeling focuses on water and solute exchanges through the matrix and conduit domain interface. To study the accuracy of the modeling approaches, we compare the CCPF and SD modeling simulations with laboratory experiment results. We also conduct a sensitivity study on the two parameters, aex and a, used in the CCPF and SD models, respectively. Based on these studies, we make the following conclusions: 1. tte continuum pipe-flow is mathematically well-posed and regularized. tte finite elementary approximations based on the mathematical model are numerically converged. 2. In comparison with the laboratory experimental results, the SD model can reproduce the experimental results of hydraulic head distribution and plume evolution. tte CCPF model overestimates the hydraulic heads along the interface between matrix and conduit, underestimates solute transport in conduits, and does not completely capture the matrix plume distribution. 3. In the CCPF model, there is a parameter, mass exchange rate, aex, to relate the water exchange between the matrix and conduit domains. The hydraulic head and solute distributions are very sensitive to the variation of the parameter. This parameter is very difficult to obtain. It is normally obtained through curve fitting using field experiments. In the SD model, there is also a parameter a in the Beavers-Joseph interface condition. However, the calculated hydraulic head and solute distributions in the system are not sensitive to the parameter variation, and so a is a dummy parameter. Therefore, the SD model does not require a parameter to describe the mass exchange between the two domains. The Stokes-Darcy modeling method can provide an independent calculation method to estimate the aex used in CCPF. ACKNOWLEDGMENTS ttis work is supported by the CMG program of the National Science Foundation under grant numbers DMS-0620035. The author wish to thank Dr. Steve Worthing- ton, Dr. Joe Donoghue, two anonymous reviewers and the special editors, Dr. Nico Goldscheider and Dr. Nataša Ravbar, for their constructive comments and revisions. REFERENCES Barenblatt, G., Zheltov, Iu. & I. Kochina, 1960: Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks.- Journal of Applied Mathematics and Mechanics, 24, 5, 1286-1303. Bauer, S., Liedl, R. & M. Sauter, 2000: Modelling of karst development considering conduit-matrix exchange flow.- In: Stauffer, F. et al. (eds.) Calibration and reliability in groundwater modelling: coping with uncertainty, Proceedings of the ModelCARE'99 Conference, 20'h-23'h September, 1999. Wallingford, Oxfordshire, UK. IAHS Publ., 265, 10-15. Bauer, S., Liedl, R. & M. Sauter, 2003: Modeling of karst aquifer genesis: Influence of exchange flow.- Water Resources Research, 39, 10, 1285, doi:10.1029/ 2003WR002218. Beavers, G. & D. Joseph, 1967: Boundary conditions at a naturally permeable wall. - Journal of Fluid Mechanics, 30, 197-207. Birk, S., Liedl, R., Sauter, M. & G. Teutsch, 2003: Hydraulic boundary conditions as a controlling factor in karst genesis.- Water Resources Research, 39, 1, 1004, doi:10.1029/2002WR001308. Bobok, E., 1993: Fluid Mechanics for Petroleum Engineers.- Elsevier Publication, pp. 412, New York. Cao, Y., Wang, H. & X. Xie, 1988: Dual-media flow models of karst areas and their application in north China.- In: Zhang, H. (ed.) Karst Hydrogeology and Karst Environment Protection, 21s' Congress of the International Association of Hydrogeologists, 10*^-15«^ October, Guilin, China. IAHS publication, 704, Oxfordshire. Chen, Y. & J. Bian, 1988: The media and movement of karst water.- In: Zhang, H. (ed.) Karst Hydrogeology and Karst Environment Protection, 21s' Congress of the International Association of Hydrogeologists, 10'h-15'h October, Guilin, China. IAHS publication, 555-564, Oxfordshire. Clemens, T., Hückinghaus, D., Sauter, M., Liedl, R. & G. Teutsch, 1996: A combined continuum and discrete network reactive transport model for the simulation of karst development.- In: Kovar, K. & P. van der Heijde (eds.) Calibration and Reliability in Groundwater Modelling, Proceedings of the ModelCARE 96, Conference, 24'h-26'h September, Golden, Colorado, USA. IAHS Publication, 237, 309-318, Oxfordshire. Faulkner, J., Hu, B. X., Kish, S. & F. Hua, 2009: Laboratory analog and Numerical study of groundwater flow and solute transport in a karst aquifer with conduit and matrix domains.- Journal of Contamination Hydrology, 110, 34-44. Ford, D. C., 2003: Perspectives in karst hydrology and cavern genesis.- Speleogenesis and Evolution of Karst Aquifers, 1, 1, 1-13. Ford, D. C. & P. W. Williams, 1989: Karst Geomorphology and Hydrology.- Academic Division of Unwin Hyd-man Ltd, pp. 601, London. Gerke, H. & M. van Genuchten, 1993a: A dual-porosity model for simulating the preferential movement of water and solutes in structured porous media.- Water Resources Research, 29, 305-319. Gerke, H. & M. van Genuchten, 1993b: Evaluation of a first-order water transfer term for variably saturated dual-porosity models.- Water Resources Research, 29, 1225-1238. Harbaugh, A., 2005: MODFLOW-2005, the U.S. geological survey modular ground-water model - The ground-water flow process.- Technical Report, USGS Techniques and Methods, no. 6-A16. Katz, B., Catches, J., Bullen, T. & R. Michel, 1998: Changes in the isotopic and chemical composition of ground water resulting from a recharge pulse from a sinking stream.- Journal of Hydrology, 211, 178207. Kiraly, L., 1998: Modeling karst aquifers by the combined discrete channel and continuum approach.- Bulletin d'Hydrogeologie, 16, 77- 98. Li, G., Loper, D. & R. Kung, 2008: Contaminant sequestration in karstic aquifers: Experiments and quantification.- Water Resources Research, 44, doi:10.1029/2006WR005797. Liedl, R., Sauter, M., Hückinghaus, D., Clemens, T. & G. Teutsch, 2003: Simulation of the development of karst aquifers using a coupled continuum pipe flow model.- Water Resources Research, 39, 3, 1057, doi:10.1029/2001WR001206. MacQuarrie, K. & E. Sudicky, 1996: On the incorporation of drains into three-dimensional variably saturated groundwater flow models.- Water Resources Research, 32, 2, 447-482. Martin, J. & R. Dean, 1999: Temperature as a natural tracer of short residence times for ground water in karst aquifers.- In: Palmer A. et al. (eds.) Karst modeling: proceedings of the symposium, 24'h-27'h Febuary 1999, Charles Town, West Virginia, USA. Karst Waters Institute Special publication, 5, 236242. Martin, J. & R. Dean, 2001: Exchange of matrix and conduit water with examples from the Floridan aquifer.- Chemical Geology, 179, 145-165. Paulson, R., Chase E., Roberts, R. & D. Moody, 1990: National Water Summary 1988-89: hydrologic Events, floods and droughts.- Technical Report, USGS Water Supply Paper Report, no. 2375. Saffman, P., 1971: On the boundary condition at the interface of a porous medium.- Stud. Appl. Math., 1, 77-84. Sauter, M., 1992: Quantification and forecasting of regional groundwater flow and transport in a karst aquifer (Gallusquelle, Malm SW Germany).- PhD thesis. Univ. of Tübingen, Germany, Tübinger Geowiss. Arb. C13, pp. 150. Shoemaker, W. B., Kuniansky, E. L., Birk, S., Bauer, S. & E. D. Swain, 2008: Documentation of a Conduit Flow Process (CFP) for MODFLOW-2005.- U.S. Geological Survey. U.S. Geological Survey Techniques and Methods, 6-A24. Taylor, C. & E. Greene, 2001: Quantitative Approaches in Characterizing Karst Aquifers.- U.S. Geological Survey Karst Interest Group Proceedings, US Geological Survey Water Resources Investigations Report 01-4011, pp. 164-166. Teutsch, G., 1989: Groundwater models in karstified terraines-two practical examples from the Swabi-an Alb, S. Germany.- In: 4"' Conference on Solving Groundwater Problems with Models, 7'h-9'h February 1989, Indianapolis, USA. International Ground Water Modeling Center, 929-953, Indianapolis, Indiana. Worthington, S., 2003: A comprehensive strategy for understanding flow in a carbonate aquifer.- Speleo-genesis and Evolution of Karst Aquifers, 1, 1, 1-8. Worthington, S., 2007: Groundwater residence times in unconfined carbonate aquifers.- Journal of Cave and Karst Studies, 69 1, 94-102.