# Turbulent and Laminar Flow in Karst Conduits Under Unsteady Flow Conditions: Interpretation of Pumping Tests by Discrete Conduit-Continuum Modeling

Markus Giese, T. Reimann, Vincent Bailly-Comte, J.-C. Marechal, M. Sauter, T. Geyer

<figure><img src="image_1.png" /></figure>

## To cite this version:

Markus Giese, T. Reimann, Vincent Bailly-Comte, J.-C. Marechal, M. Sauter, et al., Turbulent and Laminar Flow in Karst Conduits Under Unsteady Flow Conditions: Interpretation of Pumping Tests by Discrete Conduit-Continuum Modeling. Water Resources Research, 2018, 54 (3), pp.1918-1933. 10.1002/2017WR020658. [hal-02860575](https://HAL Id: hal-02860575)

Submitted on 31 Aug 2020

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Turbulent and Laminar Flow in Karst Conduits Under Unsteady Flow Conditions: Interpretation of Pumping Tests by Discrete Conduit-Continuum Modeling

M. Giese¹,²,³ ID, T. Reimann⁴ ID, V. Bailly-Comte³ ID, M. Sauter¹ ID, and T. Geyer⁵

¹Geoscientific Centre, University of Göttingen, Göttingen, Germany, ²Now at Department of Earth Sciences, University of Gothenburg, Gothenburg, Sweden, ³NRE, BRGM, University of Montpellier, Montpellier, France, ⁴Institute for Groundwater Management, TU Dresden, Dresden, Germany, ⁵Landesamt für Geologie, Rohstoffe und Bergbau im Regierungspräsidium Freiburg, Freiburg, Germany

**Abstract** Due to the duality in terms of (1) the groundwater flow field and (2) the discharge conditions, flow patterns of karst aquifer systems are complex. Estimated aquifer parameters may differ by several orders of magnitude from local (borehole) to regional (catchment) scale because of the large contrast in hydraulic parameters between matrix and conduit, their heterogeneity and anisotropy. One approach to deal with the scale effect problem in the estimation of hydraulic parameters of karst aquifers is the application of large-scale experiments such as long-term high-abstraction conduit pumping tests, stimulating measurable groundwater drawdown in both, the karst conduit system as well as the fractured matrix. The numerical discrete conduit-continuum modeling approach MODFLOW-2005 Conduit Flow Process Mode 1 (CFPM1) is employed to simulate laminar and nonlaminar conduit flow, induced by large-scale experiments, in combination with Darcian matrix flow. Effects of large-scale experiments were simulated for idealized settings. Subsequently, diagnostic plots and analyses of different fluxes are applied to interpret differences in the simulated conduit drawdown and general flow patterns. The main focus is set on the question to which extent different conduit flow regimes will affect the drawdown in conduit and matrix depending on the hydraulic properties of the conduit system, i.e., conduit diameter and relative roughness. In this context, CFPM1 is applied to investigate the importance of considering turbulent conditions for the simulation of karst conduit flow. This work quantifies the relative error that results from assuming laminar conduit flow for the interpretation of a synthetic large-scale pumping test in karst.

## 1. Introduction

Pumping tests and the interpretation of pressure or drawdown curves are frequently applied and essential tools for solving petroleum engineering or hydrogeological problems. In general, pumping tests are used to assess the hydraulic characteristics of aquifer systems (Gringarten, 1982; Kruseman & de Ridder, 1991). Besides the prediction of the overall system behavior, a wide range of specialized interpretation methods are available to quantify well and aquifer specific parameters (e.g., Bourdet et al., 1983).

Especially during the last few decades, fractured rock aquifers became the focus of attention as potential groundwater resources (e.g., Guiheneuf et al., 2014; Leray et al., 2013; Nastev et al., 2004) and waste repositories (e.g., Follin et al., 2014; Joyce et al., 2014; MacQuarrie & Mayer, 2005; Tsang et al., 2015). A wide variety of analytical techniques were developed to characterize fractured rock aquifers by pumping test evaluation (see also Agarwal et al., 1970; Bourdet, 2001; Gringarten, 1982). Karstified aquifers, as a specific group of fractured aquifer systems, display considerable complexity due to the large contrast in hydraulic parameters within the coupled conduit-matrix system. Even among karst aquifers, hydraulic parameters of conduit, fissures, and the porous matrix may vary by several orders of magnitude from local to catchment scale (e.g., Király, 2002; Sauter, 1992). As one example, the tabulation of traditionally analyzed aquifer tests and specific capacity tests for karstified Floridan Aquifer systems (mixed areas of large mature conduits and areas of preferential flow layers—large horizontal dissolution features) indicate a range of 6 orders of magnitude (Kuniansky & Bellino, 2016). The highly permeable karst conduits, draining the fissured rock matrix, are the most important hydraulic features adding a fast flow component to the groundwater discharge. Depending
on the degree of karstification, a variety of karst conduit systems with different hydraulic properties prevail. Quinlan and Ewers (1985) divide karst systems into three categories: (1) diffuse karst systems with slightly developed karstic features and mainly laminar matrix flow, (2) mixed flow karst systems with laminar matrix flow and nonlinear (turbulent) flow in solution-enlarged structures, and (3) conduit flow karst systems, where predominantly turbulent flow conditions prevail in the mature conduit systems. Especially, the conduit diameter can range from a few centimeters in initially enlarged fractures (e.g., White, 2002) up to predominant conduit structures with a large diameter of quasi-infinite hydraulic conductivity (e.g., Maréchal et al., 2008). Except for carbonate aquifers with a rather vuggy porosity as preferential flow path, randomly located individual boreholes are likely to miss these highly permeable features and can, therefore, only represent the hydraulic parameters in the vicinity of the wellbore, i.e., that of the fissured/fractured matrix (Sauter, 1992; White, 2002; Worthington, 2009). Maréchal et al. (2008) showed that long-term groundwater abstraction, with defined abstraction rate, directly from the conduit system can be a useful tool to characterize karst aquifers on catchment scale, e.g., to derive general flow patterns from the analysis of diagnostic plots (Giese et al., 2017).

Plotting pressure or drawdown curves and their derivatives with respect to time on a log-log graph (diagnostic plots) is a useful tool to obtain qualitative (identification of dominant flow regimes at specific times) and quantitative (estimation of hydraulic parameters) information about an aquifer system (Gringarten, 1982). A frequently applied method is type-curve matching where pressure or drawdown curves are compared to a set of analytical model solutions and the best fit is chosen to assess the hydraulic parameters (Bourdet, 2001; Gringarten, 1987). One drawback of the type-curve matching technique is the ambiguity in interpretation (e.g., Kruseman & de Ridder, 1991; Renard, 2005) complicating the selection of the "correct" theoretical model, most likely to produce erroneous hydraulic parameters (Gringarten, 1982; Kruseman & de Ridder, 1991). All theoretical models have in common that they are based on ideal conditions that do not represent natural aquifer conditions (Kruseman & de Ridder, 1991). Especially, the assumption of laminar flow restricts the applicability of these solutions to karstic aquifers due to the limitations of Darcy's law in quantifying turbulent flow conditions. Darcy's law is only valid for low flow velocities with small hydraulic gradients or media with narrow openings, requirements which are not achieved for water abstraction in the vicinity of or directly from the conduit system (Kruseman & de Ridder, 1991). According to Galvão et al. (2016), pumping test analyses in karst systems based on laminar flow, e.g., Theis type-curve matching, underestimate the hydrological characteristics of the conduit system.

The dimensionless Reynolds number ($Re$), representing the ratio of inertial forces versus viscous forces, is used to indicate the actual state of flow. By increasing the flow velocity, the force of inertia also increases. Flow becomes turbulent when the inertial forces overcome the viscous forces. The transition from laminar to turbulent flow is defined by the critical Reynolds number ($N_{Re}$), a guiding value depending on the hydraulic properties of the fluid and flow media, e.g., smoothness of the grains or pore walls, pore diameter, and tortuosity of the connected pore space (Shoemaker et al., 2008a). In porous media, flow is fully turbulent at Reynolds numbers above 100 and Darcy's law is only applicable for Reynolds numbers below 10 (Bear, 1972). Pipe flow is usually considered to be gradually turbulent if the Reynolds numbers exceeds 2,000 (e.g., Dreybrodt, 1988). Conduits, as preferential flow paths in karst aquifers, are considered as pipe structures and, therefore, the onset of turbulent flow is considered for Reynold numbers of approximately $N_{Re} = 500$ (White, 2002) to $N_{Re} = 2000$ (Shoemaker et al., 2008b). Even though the break point between laminar and nonlaminar flow has to be determined for every single karst system those hydrodynamic thresholds (i.e., $N_{Re}$) are also frequently applied in predefined ranges for the numerical computation of flow in idealized representations (e.g., Reimann et al., 2011) or on a catchment scale (e.g., Halihan et al., 2000; Shoemaker et al., 2008a) of karst systems. Those predefined values normally represent break points of artificial pipe systems with a smoother wall roughness and straight courses.

The conceptual laminar pipe flow model assumes a parabolic cross sectional velocity distribution inside a circular pipe. Averaging the velocity results in the Hagen-Poiseuille equation, which is frequently used for laminar pipe and/or horizontal well bore flow (e.g., Dikken, 1990; Shoemaker et al., 2008b). For nonlaminar flow, the cross sectional velocity is rather uniformly distributed which is caused by pipe roughness. This flow condition can be described by the Colebrook-White equation. The ratio of roughness height to conduit diameter, referred to as relative roughness, can be high for karst conduits, e.g., 0.25 measured by Jeannin (2001) or even higher as concluded by Atkinson (1977). The increased mean roughness height, a
nongeometrical parameter which may be caused by debris load or collapses along the flow path result in the onset of turbulence at lower Reynolds numbers and, therefore, in an increase of energy losses along the conduit flow path. Consequently, turbulent flow increases the hydraulic gradient for a defined flow rate (Dikken, 1990). Under these conditions, specific discharge is no longer linear to the head gradient (Bourdet, 2001). Those pressure differences along the flow path can also be characterized as finite conductivity, a description both used for fracture flow (e.g., Cinco-Ley et al., 1978) and flow inside horizontal wells (e.g., Dikken, 1990). The restriction in flow also changes the entire flow pattern in terms of increasing matrix inflow in the vicinity of the wellbore already at early times of water abstraction (Bourdet, 2001; Cinco-Ley et al., 1978). During this period, the drawdown signal is a superposition of linear flow in the highly conductive feature and radial flow in the matrix resulting in bilinear flow conditions which is represented by a quarter-slope line in the diagnostic plots (Bourdet, 2001; Cinco-Ley & Samaniego-V, 1981).

The Reynolds number, describing the state of flow, and the relative roughness define the geometric friction factor ($f$). For natural karst conduit systems the friction factor is normally higher by several orders of magnitude compared to manmade pipe systems (Jeannin, 2001; Springer, 2004; Worthington & Soley, 2017). Especially in sections with debris load or collapses along the flow path, the friction factor increases (Atkinson, 1977; Peterson & Wicks, 2006). In general, the friction factor of natural karst systems can vary in a wide range between $f = 0.1$ and $f = 340$ (Jeannin, 2001). Worthington and Soley (2017) analyzed the effect of turbulent flow in different karst aquifers and concluded that turbulent flow evidently increases hydraulic head and therefore needs to be considered on catchment scale. Depending on the hydraulic properties of the karst aquifer, conduit flow is already predominantly turbulent under base flow conditions (e.g., Halihan et al., 2000; White, 2002; Worthington & Soley, 2017). Consequently, model approaches for the interpretation of general flow conditions in a karst aquifer considering laminar as well as turbulent flow are capable of better representing the flow physics and, therefore, have major advantages compared to solutions solely based on Darcian flow.

The dual flow concept of karst aquifers has been incorporated into different distributive numerical groundwater flow models (Sauter et al., 2006) using different conceptual models (e.g., single continuum: Doummar et al., 2012; Mayaud et al., 2016; double continuum: Kordilla et al., 2012; discrete conduit-continuum models: de Rooij et al., 2013; Gallegos et al., 2013; Oehlmann et al., 2015; Saller et al., 2013). One numerical discrete conduit-continuum modeling approach is CFPM1 for MODFLOW-2005 (Shoemaker et al., 2008a). CFPM1 simulates laminar and turbulent pipe flow coupled through linear head-dependent water transfer with a laminar flow matrix continuum. This allows the hydraulic simulation of complex karst aquifer systems. Integrating further hydraulic features such as fast-responding drainable storage (conduit-associated drainable storage [CADS]), CFPM1 is capable of representing the characteristic drawdown of large-scale pumping tests in karst aquifers (Reimann et al., 2014).

Depending on the degree of karstification, different types of conceptual models are being applied. Kovács et al. (2005) distinguish between matrix restricted flow for mature and conduit influenced flow models for karst systems with initially enlarged flow features, depending on the general baseflow pattern of the aquifer system. Both conceptual models can be separated by a threshold depending on parameters describing the degree of heterogeneity (i.e., hydraulic and geometric properties). Currently, literature provides only a few case studies examining the influence of the hydraulic conduit properties on the flow behavior on catchment scale. Although the conduits can (and likely do) dominate the flow on catchment scale the published case studies normally do not focus on the hydraulic properties of conduits especially on relative roughness. Peterson and Wicks (2006) quoted that slight changes of the roughness significantly affect the simulated spring discharge as well as solute transport. Regardless, many case studies generally assume constant hydraulic roughness along conduits (Saller et al., 2013) or argue an insensitivity with respect to hydraulic head, spring discharge and residence time for a well-developed karst aquifer (Gallegos et al., 2013). Oehlmann et al. (2015) applied a variable roughness coefficient linearly coupled to the conduit diameter, with a positive roughness slope toward the spring. Because of the absence of a single large diameter conduit at the moderately karstified study area the authors interpret the calibrated hydraulic conduit properties as a lumped value of a conduit bundle.

The consideration of turbulent flow in conduits interacting with a discrete matrix continuum has already been addressed in several studies (e.g., Reimann et al., 2011; Shoemaker et al., 2008a). In this study, drawdown differences as well as flow pattern changes resulting from the application of laminar and turbulent flow of numerically idealized karst conduits are presented and discussed with respect to the maturity of karst aquifer systems. Analytical equations for laminar and turbulent conduit flow are employed to explain
the respective differences in head loss without any knowledge about the systems' critical Reynolds number. Therefore, the relative error caused by the application of laminar flow equations for different karst system representations is quantified. Thanks to the applied discrete conduit-continuum approach, simulations of matrix drawdown restricted by either the hydraulic conduit or the hydraulic matrix parameter can be compared. Therefore, this paper provides a guideline for the necessity of the application of turbulent flow equations in the interpretation of general flow pattern in different categories of karst systems. Another focus is set on the general problem of the application of the Reynolds number on karstic flow at regional scale. Commonly, information about the mean flow velocity of karst conduits is available from tracer experiments. In contrast, the conduit system dimension (i.e., volume, diameter, and length), which is needed for the calculation of the Reynolds number, is difficult to determine on a regional scale. Normally numerical models demand critical Reynolds numbers to discriminate between laminar and turbulent flow regimes. Due to the complexity of the hydraulic properties of karst conduits the application of such values on regional scale cannot be proven. Therefore, the study also presents the differences in head losses as a function of the hydraulic conduit properties for a defined flow velocity (i.e., drawdown signal).

## 2. Methods and Applied Tools

### 2.1. Numerical Solution: Discrete Conduit-Continuum Model CFPM1

The applied numerical method to simulate laminar Darcian flow in the fractured porous matrix continuum is based on the block-centered finite difference groundwater flow model MODFLOW-2005 (Harbaugh, 2005) according to

$$ \frac{\partial}{\partial x} \left( K_{xx} \frac{\partial h_m}{\partial x} \right) + \frac{\partial}{\partial y} \left( K_{yy} \frac{\partial h_m}{\partial y} \right) + \frac{\partial}{\partial z} \left( K_{zz} \frac{\partial h_m}{\partial z} \right) \pm W = S_s \left( \frac{\partial h_m}{\partial t} \right) \quad (1) $$

with K the hydraulic conductivity in x, y, and z direction [$L \ T^{-1}$], $h_m$ the head in the matrix cell [$L$], W the volumetric flux per unit volume [$T^{-1}$], $S_s$ the specific storage [$L^{-1}$], and time t [$T$].

The discrete conduit-continuum model CFPM1 (Shoemaker et al., 2008a) couples a discrete conduit network consisting of nodes connected by cylindrical pipes to the MODFLOW-2005 continuum. Head loss along the pipe $\Delta h_c$ [$L$] is computed by the Darcy-Weisbach equation:

$$ \Delta h_c = f \frac{\Delta l_p v^2}{2 g d_p} \quad (2) $$

with f the friction factor, $\Delta l_p$ the length of pipe [$L$], v the mean flow velocity [$L \ T^{-1}$], g the gravitational acceleration [$L^{-2}$], and $d_p$ the pipe diameter [$L$]. The Hagen-Poiseuille equation describes laminar flow as (Shoemaker et al., 2008a)

$$ Q_p = - \frac{\pi d_p^4 g \Delta h_{c, lam}}{128 v \Delta l_p \tau} \quad (3) $$

with $Q_p$ the volumetric flow rate [$L^3 \ T^{-1}$], v the kinematic viscosity of water [$L^2 \ T^{-1}$], $\Delta h_{c, lam}$ the laminar head loss [$L$], and $\tau$ the tortuosity of the pipe. Turbulent flow is computed based on the Colebrook-White equation as (Shoemaker et al., 2008a)

$$ Q_p = - \sqrt{\frac{\Delta h_{c,turb} |gd_p|^5 \pi^2}{2 \Delta l_p \tau}} \log \left[ \frac{2.51v}{\sqrt{\frac{2|\Delta h_{c,turb}|gd_p^3}{\Delta l_p \tau}}} + \frac{k_c}{3.71d_p} \right] \frac{\Delta h_{c,turb}}{|\Delta h_{c,turb}|} \quad (4) $$

with $\Delta h_{c,turb}$ the turbulent head loss [$L$] and $k_c$ the mean roughness height [$L$].

Exchange flow $Q_{ex}$ [$L^3 \ T^{-1}$] between the pipe network and the matrix continuum is calculated by a linear quasi steady state exchange coefficient $\alpha_{ex}$ [$L^2 \ T^{-1}$] as (Barenblatt et al., 1960; Bauer et al., 2003; Shoemaker et al., 2008a)

$$ Q_{ex} = \alpha_{ex} (h_c - h_m) \quad (5) $$

with $h_c$ the conduit head [$L$]. Bauer et al. (2003) introduced the exchange coefficient $\alpha_{ex}$ as a lumped conductance term representing the hydraulic characteristics as well as the geometry of the interface between conduit and matrix as
$$ \alpha_{ex} = 2\pi r_w \Delta l_p \alpha K \tau \quad (6) $$

with $r_w$ the pipe radius [L] and $\alpha$ a factor which might be interpreted as inverse fissure spacing [L⁻¹].

## 2.2. Idealized Pumping Test Analyses

To avoid the superposition of different heterogeneities on the drawdown curve, the CFPM1 model setups are based on the general requirements of an idealized aquifer (Kruseman & de Ridder, 1991):

a. constant transmissivity within the model domain of “infinite areal” extent,

b. constant pumping rate, and

c. horizontal distribution of the hydraulic head prior to pumping.

Dimensionless parameters are used for the interpretation of diagnostic plots to keep the curves independent of the magnitude of the physical parameter (Bourdet, 2001). According to Bertrand and Gringarten (1978) and Spane and Wurstner (1993), the dimensionless terms for a homogeneous, isotropic, and confined aquifer are

$$ s_D = \left( \frac{2\pi T}{Q_p} \right) \Delta s \quad (7) $$

$$ C_D = \frac{C}{2(x_f^2 S_s \pi)} \quad (8) $$

$$ t_D = \frac{Tt}{x_f^2 S} \quad (9) $$

with $s_D$ the dimensionless drawdown, $s$ the drawdown [L], $T$ the matrix transmissivity [L² T⁻¹], $C_D$ the dimensionless wellbore storage, $C$ the wellbore storage constant [L²], $x_f$ the fracture half-length [L], $t_D$ the dimensionless time, and $S$ the matrix storativity.

Frequently used tools for pumping test analysis are diagnostic plots presenting drawdown and additional drawdown derivative displayed on a log-log graph. The dimensionless drawdown derivative $s_D'$ with respect to the natural logarithm of dimensionless time $t_D$ divided by the dimensionless wellbore storage $C_D$ is given by Bourdet et al. (1983) as

$$ s_D' = \frac{\partial s_D}{\partial \ln(t_D/C_D)} \quad (10) $$

The following analysis is focused on general flow pattern differences between laminar and turbulent conduit flow for an idealized single straight conduit and can be separated into two different parts. The first part aims at the interpretation of quantitative differences in head loss between laminar and turbulent flow equation for different conduit parameters. Therefore, the analytical flow equations introduced in section 2 are applied. This analysis assumes an isolated conduit without exchange flow with the matrix. The second part of the analysis is focused on the general flow pattern inside the conduit and the influences on matrix drawdown. This part evaluates the results of the discrete conduit-continuum model CFPM1, which also considers exchange flow with the matrix. Diagnostic plots are used to explain flow pattern differences.

## 3. Results

Due to the defined hydraulic signal introduced by the water abstraction from the conduit and the preexisting nature of the conduit system, the head loss along the conduit is the only variable of the laminar and turbulent flow equation (see equations (3) and (4)). According to equation (2), the head loss is a function of the mean flow velocity, the conduit diameter, and the friction factor. The friction factor includes information about the pipe roughness and the actual flow conditions indicated by the Reynolds number. According to the Colebrook-White equation (equation (4)), head loss is a function of total flow, pipe diameter, mean roughness height, and tortuosity. Each of these parameters will be investigated regarding the influence on the drawdown behavior of a single conduit without knowledge of the systems' critical Reynolds number.

### 3.1. Analytical Head Loss Differences

Head loss along the conduit with defined pumping rate can be calculated by the analytical Hagen-Poiseuille (laminar flow) and Colebrook-White (turbulent flow) equations. Here the differences are
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Influence of different conduit characteristics on the conduit flow processes: (a) analytical head loss differences between laminar and turbulent for an isolated conduit of 1 m length and a flow rate $Q_p = 0.25 \text{ m}^3 \text{ s}^{-1}$, (b) ratio of laminar and turbulent head loss for a conduit diameter of $d_p = 0.5 \text{ m}$ and different roughness heights, and (c) ratio of laminar and turbulent head loss related to tortuosity changes for a conduit diameter of $d_p = 0.5 \text{ m}$ with mean roughness heights of $k_c = 0.01 \text{ m}$ and $k_c = 0.25 \text{ m}$.</figcaption></figure>

calculated for a conduit length of $\Delta l = 1 \text{ m}$. A defined signal of $Q_p = 0.25 \text{ m}^3 \text{ s}^{-1}$ close to the pumping well (even distribution of an abstraction rate $Q_{abs} = 0.5 \text{ m}^3 \text{ s}^{-1}$ in both directions) is used as flow rate according to equations (3) and (4).

Figure 1a presents the head loss differences calculated by the analytical equations for laminar and turbulent flow pattern for two different conduit diameters ($d_p = 0.5 \text{ m}$ and $d_p = 2.5 \text{ m}$) plotted against flow velocity. Both conduits are artificially smooth with a mean roughness height of $k_c = 1 \times 10^{-4} \text{ m}$. Due to different cross-sectional areas, flow velocities inside the conduit differ depending on the conduit diameter. According to the basic flow velocity equation $v = Q_p/A$ with $Q_p = 0.25 \text{ m}^3 \text{ s}^{-1}$ and $A = \pi(d_p^2/4)$ the mean flow velocities for the applied conduit diameter of $d_p = 0.5 \text{ m}$ and $d_p = 2.5 \text{ m}$ are $v_{0.5} = 1.27 \text{ m s}^{-1}$ and $v_{2.5} = 0.051 \text{ m s}^{-1}$, respectively. Therefore, also the calculated Reynolds numbers differ. The Reynolds number for the conduit with a conduit diameter of $d_p = 0.5 \text{ m}$ is $Re = 5.6 \times 10^5$, whereas the Reynolds number for the larger conduit diameter can be calculated as $Re = 1.1 \times 10^5$.

Assuming a constant flow rate, the head loss differences between laminar and turbulent flow for the conduit diameter of $d_p = 2.5 \text{ m}$ are lower compared to a small diameter. The example in Figure 1a shows a difference of 4 orders of magnitude. The head loss difference along the conduit segment with a diameter $d_p = 0.5 \text{ m}$ is $h_{c,turb} - h_{c,lam} = 2.5 \times 10^{-3} \text{ m}$ whereas a diameter $d_p = 2.5 \text{ m}$ accounts for head loss
differences of $h_{c,turb} - h_{c,lam} = 5.5 \times 10^{-7}$ m. Even for the same velocity and therefore a higher abstraction rate ($Q_p = 6.2 \text{ m}^3 \text{s}^{-1}$) the head loss differences are lower by 1 order of magnitude. The same head loss differences between the two types of conduits can be derived by an increase in abstraction rate in the conduit with a diameter $d_p = 2.5$ m. To generate a similar head loss differences of $h_{c,turb} - h_{c,lam} = 2.5 \times 10^{-3}$ m the flow rate has to be increased to $Q_p = 16.7 \text{ m}^3 \text{s}^{-1}$ ($v_{2.5} = 3.4 \text{ m s}^{-1}$) with a Reynolds number of $Re = 7.4 \times 10^6$ (Figure 1a).

Figure 1b analyzes the effect of the mean roughness height on the head loss differences for a conduit diameter of $d_p = 0.5$ m. An increase in mean roughness height decreases the effective flow cross-sectional area or, to mention the effect, increases the probability of eddies and cross flow. As a result, the head loss difference along the flow paths increases. Figure 1b shows the ratio of laminar and turbulent head loss as a function of the Reynolds number for different roughness heights. Already for small Reynolds numbers head loss differences up to 20% are observed. The significance of head loss differences rises with increasing mean roughness heights. For Reynolds numbers similar to those in the above example ($Re \approx 5 \times 10^5$), all curves show significant head loss differences.

The tortuosity is a factor accounting the effective conduit length. Applying a constant flow rate, an increase in conduit length needs to be counterbalanced by an increase in head gradient. Figure 1c presents the head loss ratio for two different mean roughness heights of Figure 1b ($k_c = 0.01$ m and $k_c = 0.25$ m) for three different tortuosity values: (a) $\tau = 1$, (b) $\tau = 1.5$, and (c) $\tau = 2$. The three curves representing different tortuosity values for defined mean roughness heights are virtually congruent. Therefore, the impact of tortuosity changes on head loss is insignificant within a reasonable range of values.

## 3.2. Numerical Flow Pattern Differences

The quasi-infinite aquifer is considered by the large horizontal extent of the discrete conduit-continuum model domain to avoid the effects of the no-flow (Neumann) boundary conditions in the matrix continuum on the general flow pattern. The matrix continuum with an extent of 113,000 m × 113,000 m has a bottom elevation of 0 m with an aquifer thickness of $b = 250$ m. The area around the pumping well, located in the center of the domain, is discretized with 1 m cell length perpendicular to the conduit. The spatial discretization increases stepwise up to a cell length of 100 m. The spatial discretization along the conduit is set to a cell length of 20 m. The fissured matrix is considered as a confined layer. The hydraulic parameters of the fissured matrix are $K_m = 1 \times 10^{-4}$ m s$^{-1}$ and $S_s = 1 \times 10^{-4}$ m$^{-1}$. The well has a constant pumping rate of $Q_{abs} = 0.5 \text{ m}^3 \text{s}^{-1}$, implemented as (negative) Neumann boundary condition (CFPM1 - CRCH Package), for a duration of $T_{total} = 6 \times 10^7$ s. Groundwater temperature is set to 15°C.

All of the following setups consider a single conduit with a length of $L_p = 3,000$ m subdivided into 150 sections each with a section length of 20 m (total number of nodes: 151). The conduit nodes are located in the center of the matrix cell. No wellbore/conduit storage is assumed. The wellbore is centrally arranged at node 76. The uniform matrix-conduit exchange is considered with $\alpha K = 0.01 \text{ s}^{-1}$ (equation (6)). The roughness height is constant along the conduit with $k_c = 1 \times 10^{-4}$ m, representing a very smooth conduit, and the tortuosity for all conduit segments is set to $\tau = 1$. The applied critical Reynolds numbers secure the calculation according to either the Hagen-Poiseuille (laminar flow) or Colebrook-White (turbulent flow) equation without any transition during the model time $T_{Total}$. The interpretation of drawdown and flow processes (Figures 2 and 3) is performed for one of the conduit end nodes 1 (N1) and node 76 (N76) containing the pumping well.

### 3.2.1. Influence on Conduit Flow Pattern: Conduit Diameter

Figures 2a and 2b present diagnostic plots of two different parameterizations according to the above mentioned conduit diameters. Differences between turbulent and laminar flow regimes become apparent for a conduit diameter of $d_p = 0.5$ m. The drawdown curve derived from the turbulent flow equation (Figure 2a) shows higher conduit drawdown during all periods of pumping. Especially in the beginning of groundwater abstraction, drawdown differences are significant and the diagnostic plots do not show the typical shape of linear flow conditions. The diagnostic plot shows a quarter unit slope and therefore bilinear flow conditions can be assumed. This flow behavior is a clear evidence of reduced conduit conductivity (finite conductivity conduit). This finding is also supported by the head difference along the conduit (Figure 2g).

Figure 2c shows the conduit flow during the pumping test for the small conduit diameter. Due to the additional head losses caused by turbulent flow, flow toward the pumping well is reduced at the beginning of
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Influence of different conduit diameters on the conduit flow processes; (left) for a conduit diameter of $d_p$ = 0.5 m and (right) for a conduit diameter of $d_p$ = 2.5 m: (a, b) diagnostic plots; (c, d) conduit flow toward the pumping well; (e, f) exchange flow at the pumping well (N76) and the tip of the conduit (N1); (g, h) hydraulic head difference between pumping well (N76) and conduit tip (N1).</figcaption></figure>

pumping. As a consequence, the exchange flow in the vicinity of the pumping well is high but constantly declining with increasing conduit flow (Figure 2e). At the tips of the conduit the drawdown increases only slightly resulting in an insignificant exchange flow during this early stage of pumping (Figure 2e). Hydraulic head differences between the pumping well and the tips of the conduit create a continuously rising hydraulic gradient (Figure 2g) which increases the conduit flow toward the pumping well (Figure 2c). As a consequence, the exchange flow at the pumping well drops and the increased head difference between matrix and conduit increases at the same time the exchange flow rate between the matrix and the conduit tips (N1, Figure 2e). According to the general definition of linear flow (dominating inflow at the conduit end) the point of intersection roughly indicates a flow pattern change from bilinear to linear flow. The linear flow period lasts only for a short period before all fluxes and also the hydraulic gradient approach steady state condition (Figures 2c, 2e, and 2g). During "steady state" radial flow $\Delta h_{1,500 \text{ m}}$, the head difference between pumping well and the tip of the conduit, equals 1.4 m which amounts to a total head loss of $9.3 \times 10^{-4} \text{ m}$ per conduit meter.

Applying the laminar flow equation on the small diameter conduit ($d_p$ = 0.5 m) results in a lower head loss along the conduit. Therefore, the conduit flow toward the pumping well nearly reaches the flow maximum
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Influence of roughness height on the conduit flow processes: comparison for turbulent flow conditions for a mean roughness height of $k_c$ = 0.01 m (black curves) and $k_c$ = 0.25 m (red curves): (a) diagnostic plots; (b) conduit flow toward the pumping well (N76); (c) exchange flow at the pumping well and the tip of the conduit; (d) hydraulic head difference between pumping well (N76) and conduit tip (N1).</figcaption></figure>

at the start of pumping (Figure 2c). The exchange flow rate in the vicinity of the pumping well is lower compared to that for turbulent flow (Figure 2e). At the beginning of pumping the exchange flow rate is nearly constant along the conduit. Hence, the point of intersection between the exchange flow curves, indicating the start of linear flow, is reached early in time. During the period of radial flow condition, flow patterns are comparable to those for turbulent flow conditions. At the end of the pumping period the head difference between pumping well and the tip of the conduit is $\Delta h_{1,500 \text{ m}} = 1.8 \times 10^{-2} \text{ m}$. The head loss difference along the conduit between turbulent and laminar flow conditions is $\Delta h_{c,turb} - \Delta h_{c,lam} = 9.2 \times 10^{-4} \text{ m}$ for the given parameterization. This value approximately equals that of the analytical solution considering conduit flow without matrix exchange. The exchange flow with the matrix accounts for the minor differences between the analytical and numerical result.
The laminar as well as the turbulent flow behavior of the conduit with a conduit diameter $d_p = 2.5$ m is comparable to that of the laminar flow equation for the small diameter conduit ($d_p = 0.5$ m). At the start of pumping the exchange flow is constant along the conduit and both diagnostic plots already indicate linear flow (Figure 2b). Drawdown inside the conduit is uniform during the entire duration of pumping which results in negligible head differences along the conduit (Figure 2h). This characterization of the flow components based on drawdown curves analyses shows that the head loss difference between laminar and turbulent flow conditions is negligible. This also fits the results of the analytical flow equations (Figure 1a). The hydraulic conduit parameters do not affect the linear flow pattern. Due to the uniform matrix exchange, flow is only restricted by the hydraulic properties of the interface and the matrix.

### 3.2.2. Influence on Conduit Flow Pattern: Mean Roughness Height

The influence of roughness height changes on the general flow behavior is presented in Figure 3. Figure 3a shows a diagnostic plot for mean roughness heights of $k_c = 0.01$ m and $k_c = 0.25$ m (cf., Figure 1b) for a conduit diameter of $d_p = 0.5$ m. The other parameters are similar to those introduced in section 3.2.

Compared to the smooth conduit of section 3.2.1, a raise of the mean roughness height influences the conduit flow properties, e.g., increase the turbulent core zone. The interferences caused by the wall roughness, decrease the conductivity and thus the flux along the conduit (Figure 3b). Therefore, the effect of flow restriction is further increased, causing a higher exchange flow rate in the vicinity of the pumping well (Figure 3c) as well as increased matrix inflow near the pumping well (Figure 3d). The effects of increased mean roughness height on conduit drawdown are similar to those already explained above in Figure 2. One exception is presented by the drawdown curve scenario with a mean roughness height of $k_c = 0.25$ m. With a roughness height equal to the radius of the conduit ($k_c/d_p = 0.5$), the flow restriction along the conduit is so high that the flux along the conduit is less than half of that of the initial setup. Due to reduced conduit flow, groundwater is mainly abstracted from the matrix resulting in radial flow conditions already during early times of the pumping test (Figure 3a).

### 3.3. Influence of Conduit Drawdown on Matrix Heads

Beside the differences in conduit drawdown also differences in matrix drawdown can be observed caused by turbulent flow in small diameter conduits. Figure 4 shows the matrix drawdown of a 5,000 m × 5,000 m area around the conduit center (node 76) during (bi)linear flow, transition period and radial flow for the model parameterization described in section 3.2 ($d_p = 0.5$ m). Figure 4 presents the matrix drawdown at the time steps marked in Figure 2a.

For laminar flow condition, the head losses inside the conduit are low resulting in a more or less uniform drawdown along the conduit. As a consequence, the exchange flow from the adjacent matrix (cf., Figure 2e) as well as the matrix drawdown is nearly uniform along the conduit (Figure 4, laminar).

The head distribution shows significant spatial differences for turbulent flow conditions. These differences become apparent in conduit as well as matrix heads. The quadratic conduit head losses, caused by turbulence, result in a distinctive head gradient along the conduit (Figure 2g). Flow restrictions affect conduit flux as well as the exchange flow rate between conduit and matrix. As a consequence of the decreased flux, the drawdown at the pumping well increases (Figure 2a). The reduced conduit flux causes a high matrix flux toward the pumping well. Hence, matrix exchange flow (Figure 2e) and the resulting matrix drawdown are nonuniform along the conduit (Figure 4). Linear conduit flow superimposed by radial flow inside the adjacent rock toward the well generates bilinear flow conditions. The hydraulic parameters of the conduit clearly influence the flow pattern of the matrix and, therefore, flow can be considered as conduit influenced. Compared to laminar flow conditions, the drawdown in the vicinity of the pumping well is higher, decreasing with distance to the well toward the tips of the conduit. Because of the reduced flux the conduit drawdown at the conduit end is lower than for laminar flow conditions. This affects matrix exchange flow and hence results in a reduced matrix drawdown. The matrix drawdown differences are also presented in Figure 4. Positive head differences mark the location where the turbulent conduit flow increases the matrix drawdown and negative head differences were obtained in locations where laminar flow increases the matrix drawdown.

The general flow behavior does not change until the start of the transition period. The differences in matrix drawdown adjacent to the pumping well further increase as well as the differences at the conduit tips. For linear flow, groundwater abstraction from the matrix is still uniform along the conduit. With increasing
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Matrix head distribution along a horizontal conduit ($d_p$ = 0.5 m) during large-scale constant groundwater abstraction at the time steps marked in Figure 2. (left) Head distribution along a conduit with laminar flow; (middle) head distribution along a conduit with turbulent flow; (right) drawdown differences between the two flow conditions showing higher spatial matrix drawdown during turbulent flow in red and higher spatial matrix drawdown of laminar flow in blue.</figcaption></figure>

pumping duration and expanding cone of depression, the differences of flow pattern between laminar and turbulent flow conditions vanish, but the overall drawdown difference increases.

## 4. Discussions

According to the results of section 3, head losses and flow pattern are highly related to the hydraulic properties of the conduit. Literature provides different critical Reynolds numbers for the transition between laminar to turbulent conduit flow even the break point between laminar and turbulent flow needs to be determined by physical experiments. Figure 5 shows the analytical head loss differences as well as the calculated friction factor as a function of the mean roughness for a single conduit with a diameter $d_p$ = 0.5 m. The calculations are based on two different Reynolds numbers: $Re = 500$ (Figure 5a, $Q_p = 2.2 \times 10^{-4}$ m³ s⁻¹, $v = 1.1 \times 10^{-3}$ m s⁻¹) and $Re = 10,000$ (Figure 5b, $Q_p = 4.5 \times 10^{-3}$ m³ s⁻¹, $v = 2.3 \times 10^{-2}$ m s⁻¹).

The results presented in Figure 5a can be divided at least into two different parts. Below a relative roughness of $k_c/d_p$ = 0.01, the head loss differences are low and uniform. Also the calculated friction factor (based on equation
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Head loss difference between laminar and turbulent flow and friction factor for different hydraulic conduit parameter for a flow rate $Q_p = 0.25 \text{ m}^3 \text{ s}^{-1}$: (a) head loss difference and friction factor for $Re = 500$ ($d_p = 0.5$ m), (b) head loss difference and friction factor for $Re = 10,000$ ($d_p = 0.5$ m), (c) head loss difference for different conduit diameters and mean roughness heights, and (d) friction factor for different conduit diameters and mean roughness heights.</figcaption></figure>

(2)), divided into a laminar and turbulent portion is nearly constant. Due to the independence of the laminar head loss from the mean roughness height, the laminar friction factor is constant along the abscissa. Starting at a relative roughness of $k_c/d_p = 0.01$, the turbulent friction factor steadily increases. At a relative roughness of $k_c/d_p = 1$, the turbulent friction factor is $f = 0.75$ but the head loss difference $h_{c,turb} - h_{c,lam} = 1 \times 10^{-7} \text{ m}$ is still insignificant. In Figure 5b, the head loss differences are generally higher because of the increased flow rate. Nevertheless, the friction factor for a relative roughness of $k_c/d_p = 1$ is comparable to that for a Reynolds number of $Re = 500$. Therefore, this part can be referred to as hydraulic rough. The friction factor only depends on the relative roughness. Hence, the range between $k_c/d_p = 0.01$ and $k_c/d_p = 1$ describes the transition period between laminar and turbulent flow. For a Reynolds numbers $Re = 10,000$ significant head loss differences, able to change the flow pattern, can only be derived by friction factors higher than approximately $f = 10$. According to Figure 5b, these friction factors can only be achieved by a relative roughness beyond $k_c/d_p = 1$.
Normally, for regional scale the hydraulic properties of karst conduits are partly or totally unknown. Figure 5c presents head loss differences between laminar and turbulent flow conditions and Figure 5d calculated friction factors per conduit meter related to the conduit diameter and the mean roughness height. The mean roughness height covers the whole range from artificially smooth ($k_c$ = 0.001 m) to the respective conduit diameter ($d_p$ = 1). The conduit diameter ranges between $d_p$ = 0.5 m and $d_p$ = 2.5 m. The values presented in Figures 5c and 5d are based on equations (3) and (4) with a flow rate $Q_p$ = 0.25 m³ s⁻¹ and no knowledge of the critical Reynolds number and which equation actually applies. The calculated mean velocity ranges between $v_{0.5}$ = 1.27 m s⁻¹ and $v_{2.5}$ = 0.05 m s⁻¹ and is therefore higher than that applied in the example presented in Figures 5a and 5b.

Slightly developed karst systems, e.g., the Gallusquelle catchment (e.g., Oehlmann et al., 2015; Sauter, 1992), commonly do not have apparent large conduit structures. Based on the results of a calibrated distributed parameter flow and transport model, Oehlmann et al. (2015) concluded that the surface to volume ratio is high for the conduit network. Flow is likely to be dominated by bundles of small scale karst flow features. According to Figure 5c, the head loss differences, between laminar and turbulent flow, are highest for conduits of small diameter combined with high mean conduit roughness. Therefore, the application of a laminar instead of a turbulent flow equation for parameter estimation is likely to lead to significant errors. The analyses of large-scale pumping tests in karst systems with slightly enlarged flow features, well connected to the fissured matrix, or for systems with only a low permeability contrast would reveal a high exchange flow in the vicinity of the pumping well leading to bilinear flow (see also section 3). For conduit networks with high conduit storage, not illustrated here, a higher volume of water will be drained from the storage most probably masking the response of conduit flow during early times.

Mixed flow karst systems and mature karst systems are dominated by dissolution-enlarged conduit systems with (partly) large conduit diameter, for example, the Cent Fonts catchment (Maréchal et al., 2008). For a constant conduit diameter, the differences in head loss are insignificant, even for high mean roughness values (Figure 5c). Based on this assumption and employing linear flow equations, parameter estimations can be considered as relatively accurate in mature karst systems. This could also explain the results of Gallegos et al. (2013) in terms of the insensitivity with regard to mean roughness on subregional scale during base flow in a well-developed karst aquifer (cf., Kuniansky, 2016). The errors of the estimated parameters will be insignificant even for high pumping rates as used during the large-scale pumping test at the Cent Fonts catchment. According to the results, the approximation of quasi-infinite hydraulic conductivity by Maréchal et al. (2008) can be confirmed.

The applied setups for the numerical conduit-continuum model use an idealized parametrization. The parameter combinations are used to minimize the effects on drawdown behavior caused by processes other than the type of flow (laminar/turbulent) in the karst conduit. One of these processes is the exchange with the limestone matrix especially influencing the drawdown at the beginning of pumping. Another simplification is related to the storativity. Changing the matrix storativity has a negligible influence on the drawdown curve and does not influence the general flow pattern. Furthermore, all setups do not consider fast-responding storage in karst conduits, which mask the drawdown behavior at the beginning of pumping. The effect of fast-responding conduit storage is already described by Reimann et al. (2014) and Giese et al. (2017).

Additionally, the above stated results for turbulent flow in karst conduits are computed with the Colebrook-White equation. This equation applies to pressurized flow at moderate Reynolds numbers (transition zone of laminar and turbulent flow) and small diameter pipes with natural roughness. Increased roughness, for example, due to deposits or the natural shape of karst conduits, are not considered. Therefore, the use of the Colebrook-White equation already idealized flow conditions and the applicability on certain karst aquifer systems must be examined in detail.

# 5. Conclusions

The above analysis shows the need for the consideration of turbulent flow in karst aquifer modeling and characterization especially for those aquifer systems defined as karst systems with slightly enlarged flow features. Turbulent flow may result in restricted flow inside the conduit, also referred to as finite conductivity. The influence of turbulent flow conditions on drawdown is especially large for a high relative roughness
(small conduit radii and high mean roughness height). For those conditions, turbulent flow cannot be neglected in the simulation of flow physics. Otherwise the extent of water level drawdown will be overestimated at the beginning of pumping. Hence, assuming laminar flow conditions will result in an underestimation of the conduit dimension. The determination of information with respect to the conduit geometric and hydraulic properties, which is always a critical aspect in karst aquifer characterization, is required. Due to nonlinearity of hydraulic head and discharge using turbulent flow equations it is not possible to apply only one (dimensionless) parameter, e.g., $k_c/d_p$ describing the head losses during turbulent flow. As a consequence, the parameters describing the conduit geometry, especially the roughness, will serve as a calibration parameter for numerical models applied in slightly developed karst aquifers. For mature karst systems with well-developed tertiary porosity (conduit systems) approaches applying laminar flow equations will be sufficient. The analysis proofs for different Reynolds numbers and flow rates that the relative roughness must be high for causing significant head loss differences between laminar and turbulent flow pattern. Specifically for conduit systems with less developed connectivity to the adjacent fissured matrix, this assumption is adequate. Mature karst systems with collapses or high debris load along the preferential flow path can be a possible exception. For those systems, relative roughness and the friction factor are high (Atkinson, 1977). The other possible exceptions might be systems at high flow conditions, even though observations of karst systems with broad conduit diameter (e.g., Wakulla spring, Cent Fonts) show laminar flow behavior for Reynolds numbers clearly indicating turbulent flow.

Different studies show that the friction factor for preferential flow path in karst aquifers can be higher than $f=1$ (Jeannin, 2001). According to Figure 5d, those values can only be the result of relative roughness higher than $k_c/d_p = 1$. Springer (2004) provides one approach to separate the total head loss into three different origins: (a) skin head loss, (b) head loss as consequence of expansion, and (c) head loss caused by flow orientation. The local scale analysis of a cave reveals that the head loss caused by skin effects is comparable low to the other two. Macroscopic channel expansions and bends, in general cross-section changes, account for major head losses (Springer, 2004; Worthington & Soley, 2017). Especially for the large-scale modeling of karst systems the separation of the total head loss in different trigger is difficult. Therefore, it seems reasonable to consider the roughness as a lumped parameter reflecting roughness as well as geometrical conduit properties of the collection of preferential flow paths.

## Acknowledgments

This project was funded by the Deutsche Forschungsgemeinschaft (DFG) under grants LI 727/11-2 and GE 2173/2-2 and by the BRGM under grants PDR13D3E91 and PDR14D3E61. The authors thank Harihar Rajaram, Eve Kuniansky, Zexuan Xu, and one anonymous reviewer for constructive comments that significantly improved the manuscript. We also express our appreciation to Rudolf Liedl from the TU Dresden for his internal review. All input data for the CFPM1 are available for download: https://doi.org/10.5281/zenodo.1122635.

## References

Agarwal, R. G., Al-Hussainy, R., & Ramey, H. J. (1970). An investigation of wellbore storage and skin effect in unsteady liquid flow: I. Analytical treatment. *SPE Journal*, 10(3), 279-290.

Atkinson, T. C. (1977). Diffuse flow and conduit flow in limestone terrain in the Mendip Hills, Somerset (Great Britain). *Journal of Hydrology*, 35, 93-110.

Barenblatt, G. I., Zheltov, I. P., & Kochina, I. N. (1960). Basic concepts in the theory of seepage of homogeneous liquids in fissured rock. *Journal of Applied Mathematics and Mechanics (PMM)*, 24(5), 1286-1303.

Bauer, S., Liedl, R., & Sauter, M. (2003). Modeling of karst aquifer genesis: Influence of exchange flow. *Water Resources Research*, 39(10), 1285. https://doi.org/10.1029/2003WR002218

Bear, J. (1972). *Dynamics of fluid in porous media*. New York, NY: Dover Publications.

Bertrand, L., & Gringarten, A. C. (1978). *Détermination des caractéristiques hydrauliques des aquifères fissurés par pompage d'essai en régime transitoire—Application aux nappes de la craie* (Final Rep. 78 SGN 669 GEG) [in French], Orléans, France: BRGM.

Bourdet, D. (2001). *Well test analysis: The use of advanced interpretation models*. Amsterdam, the Netherlands: Elsevier.

Bourdet, D., Whittle, T. M., Douglas, A. A., & Pirard, Y. M. (1983). A new set of type curves simplifies well test analysis. *World Oil*, 196, 95-106.

Cinco-Ley, H., & Samaniego-V, F. (1981). Transient pressure analysis for fractured wells. *SPE Journal*, 7490, 1749-1766.

Cinco-Ley, H., Samaniego-V, F., & Dominguez-A, N. (1978). Transient pressure behaviour for a well with a finite-conductivity vertical fracture. *SPE Journal*, 6014, 253-264.

de Rooij, R., Perrochet, P., & Graham, W. (2013). From rainfall to spring discharge: Coupling conduit flow, subsurface matrix flow and surface flow in karst systems using a discrete-continuum model. *Advances in Water Resources*, 61, 21-41. https://doi.org/10.1016/j.advwatres.2013.08.009

Dikken, B. J. (1990). Pressure drop in horizontal wells and its effect on production performance. *SPE Journal*, 19824, 1426-1433. https://doi.org/10.2118/19824-PA

Doummar, J., Sauter, M., & Geyer, T. (2012). Simulation of flow processes in large scale karst systems with an integrated catchment model (Mike She)—Identification of relevant parameters influencing spring discharge. *Journal of Hydrology*, 426-427, 112-123. https://doi.org/10.1016/j.jhydrol.2012.01.021

Dreybrodt, W. (1988). *Processes in Karst systems—Physics, chemistry, and geology*. Berlin, Germany: Springer.

Follin, S., Hartley, L., Rhen, I., Jackson, P., Joyce, S., Roberts, D., et al. (2014). A methodology to constrain the parameters of a hydrogeological discrete fracture network model for sparsely fractured crystalline rock, exemplified by data from the proposed high-level nuclear waste repository site at Forsmark, Sweden. *Hydrogeology Journal*, 22(2), 313-331. https://doi.org/10.1007/s10040-013-1080-2
Gallegos, J. J., Hu, B. X., & Davis, H. (2013). Simulating flow in karst aquifers at laboratory and sub-regional scales using MODFLOW-CFP. *Hydrogeology Journal*, **21**(8), 1749–1760. https://doi.org/10.1007/s10040-013-1046-4

Galvão, P., Halihan, T., & Hirata, R. (2016). The karst permeability scale effect of Sete Lagoas, MG, Brazil. *Journal of Hydrology*, **532**, 149–162.
https://doi.org/10.1016/j.jhydrol.2015.11.026

Giese, M., Reimann, T., Liedl, R., Maréchal, J.-C., & Sauter, M. (2017). Application of the flow dimension concept for numerical drawdown data analyses in mixed flow karst systems. *Hydrogeology Journal*, **25**(3), 799–811. https://doi.org/10.1007/s10040-016-1523-7

Gringarten, A. C. (1982). Flow-test evaluation of fractured reservoirs. In T. N. Narasimhan (Ed.), *Recent trends in hydrogeology* (Spec. Pap., Vol. 189, pp. 237–262). Boulder, CO: Geological Society of America.

Gringarten, A. C. (1987). Type-curve analysis: What it can and cannot do. *Journal of Petroleum Technology*, **39**(1), 11–13.

Guiheneuf, N., Boisson, A., Bour, O., Dewandel, B., Perrin, J., Dausse, A., et al. (2014). Groundwater flows in weathered crystalline rocks: Impact of piezometric variations and depth-dependent fracture connectivity. *Journal of Hydrology*, **511**, 320–334. https://doi.org/10.1016/j.jhydrol.2014.01.061

Halihan, T., Sharp, J., & Mace, R. E. (2000). Flow in the San Antonio segment of the Edwards aquifer: Matrix, fractures, or conduits? In I. D. Sasowski & C. M. Wicks (Eds.), *Groundwater flow and contaminant transport in carbonate aquifers* (pp. 129–146). Rotterdam, the Netherlands: Balkema.

Harbaugh, A. W. (2005). *MODFLOW-2005, The U.S. Geological Survey modular ground-water model—The ground-water flow process* (Tech. Methods 6-A16). Reston, VA: U.S. Geological Survey.

Jeannin, P. Y. (2001). Modelling flow in phreatic and epiphreatic karst conduits in the Hölloch cave (Muotatal, Switzerland). *Water Resources Research*, **37**(2), 191–200.

Joyce, S., Hartley, L., Applegate, D., Hoek, J., & Jackson, P. (2014). Multi-scale groundwater flow modeling during temperate climate conditions for the safety assessment of the proposed high-level nuclear waste repository site at Forsmark, Sweden. *Hydrogeology Journal*, **22**(6), 1233–1249. https://doi.org/10.1007/s10040-014-1165-6

Király, L. (2002). *Karstification and groundwater flow*. Paper presented at the Proceedings of the conference on evolution of karst: From pre-karst to cessation (pp. 155–190), Postojna-Ljubljana.

Kordilla, J., Sauter, M., Reimann, T., & Geyer, T. (2012). Simulation of saturated and unsaturated flow in karst systems at catchment scale using a double continuum approach. *Hydrology and Earth System Sciences*, **16**, 3909–3923. https://doi.org/10.5194/hess-16-3909-2012

Kovács, A., Perrochet, P., Király, L., & Jeannin, P.-Y. (2005). A quantitative method for the characterisation of karst aquifers based on spring hydrograph analysis. *Journal of Hydrology*, **303**, 152–164. https://doi.org/10.1016/j.jhydrol.2004.08.023

Kruseman, G. P., & de Ridder, N. A. (1991). *Analysis and evaluation of pumping test data*. Wageningen, the Netherlands: Institute for Land Reclamation and Improvement.

Kuniansky, E. L. (2016). *Simulating groundwater flow in karst aquifers with distributed parameter models—Comparison of porous-equivalent media and hybrid flow approaches* (Sci. Invest. Rep. 2016–5116). Reston, VA: U.S. Geological Survey. https://doi.org/10.3133/sir201665116

Kuniansky, E. L., & Bellino, J. C. (2016). *Tabulated transmissivity and storage properties of the Floridan Aquifer system in Florida and parts of Georgia, South Carolina, and Alabama* (ver. 1.1, May 2016), U.S. Geological Survey Data Series 699.

Leray, S., de Dreuzy, J. R., Bour, O., & Bresciani, E. (2013). Numerical modeling of the productivity of vertical to shallowly dipping fractured zones in crystalline rocks. *Journal of Hydrology*, **481**, 64–75. https://doi.org/10.1016/j.jhydrol.2012.12.014

MacQuarrie, K. T. B., & Mayer, K. U. (2005). Reactive transport modeling in fractured rock: A state-of-the-science review. *Earth-Science Review*, **72**(3–4), 189–227. https://doi.org/10.1016/j.earscirev.2005.07.003

Maréchal, J.-C., Ladouche, B., Dörflinger, N., Lachassagne, P., (2008). Interpretation of pumping tests in a mixed flow karst system. *Water Resources Research*, **44**, W05401. https://doi.org/10.1029/2007WR006288

Mayaud, C., Wagner, T., Benischke, R., & Birk, S. (2016). Understanding changes in the hydrological behaviour within a karst aquifer (Lur-bach system, Austria). *Carbonates Evaporites*, **31**, 357–365. https://doi.org/10.1007/s13146-013-0172-3

Nastev, M., Savard, M. M., Lapcevic, P., Lefebvre, R., & Martel, R. (2004). Hydraulic properties and scale effects investigation in regional rock aquifers, south-western Quebec, Canada. *Hydrogeology Journal*, **12**(3), 257–269. https://doi.org/10.1007/s10040-004-0340-6

Oehlmann, S., Geyer, T., Licha, T., & Sauter, M. (2015). Reducing the ambiguity of karst aquifer models by pattern matching of flow and transport on catchment scale. *Hydrology and Earth System Sciences*, **19**, 893–912. https://doi.org/10.5194/hess-19-893-2015

Peterson, E. W., & Wicks, C. M. (2006). Assessing the importance of conduit geometry and physical parameters in karst systems using the storm water management model (SWMM). *Journal of Hydrology*, **329**, 294–305.

Quinlan, J. F., & Ewers, R. O. (1985). Ground water flow in limestone terranes: Strategy rationale and procedure for reliable, efficient monitoring of ground water quality in karst areas. In *Proceedings of the national symposium and exposition on aquifer restoration and ground water monitoring* (5th, Columbus, Ohio), pp. 197–243). Worthington, OH: National Water Well Association.

Reimann, T., Giese, M., Geyer, T., Liedl, R., Maréchal, J.-C., & Shoemaker, W. B. (2014). Representation of water abstraction from a karst conduit with numerical discrete-continuum models. *Hydrology and Earth System Sciences*, **18**, 227–241. https://doi.org/10.5194/hess-18-227-2014

Reimann, T., Rehrl, C., Shoemaker, W. B., Geyer, T., & Birk, S. (2011). The significance of turbulent flow representation in single-continuum models. *Water Resources Research*, **47**, W09503. https://doi.org/10.1029/2010WR010133

Renard, P. (2005). 160: Hydraulic of wells and well testing. In M. Anderson (Ed.), *Encyclopedia of hydrological sciences, Part 13*. Groundwater. New York, NY: John Wiley. https://doi.org/10.1002/0470848944.hsa154a

Saller, S. P., Ronayne, M. J., & Long, A. J. (2013). Comparison of a karst groundwater model with and without discrete conduit flow. *Hydrogeology Journal*, **21**(7), 1555–1566. https://doi.org/10.1007/s10040-013-1036-6

Sauter, M. (1992). Quantification and forecasting of regional groundwater flow and transport in a karst aquifer (Gallusquelle, Malm, SW Germany). *Tübinger Geowissenschaftliche Arbeiten, Part C*, **13**, 151.

Sauter, M., Kovács, A., Geyer, T., & Teutsch, G. (2006). Modellierung der Hydraulik von Karstgrundwasserleitern—Eine Übersicht [in German]. *Grundwasser*, **11**(3), 143–153.

Shoemaker, W. B., Cunningham, K. J., Kuniansky, E. L., & Dixon, J. (2008a). Effects of turbulence on hydraulic heads and parameter sensitivities in preferential groundwater flow layers. *Water Resources Research*, **44**, W03501. https://doi.org/10.1029/2007WR006601

Shoemaker, W. B., Kuniansky, E. L., Birk, S., Bauer, S., & Swain, E. D. (2008b). *Documentation of a Conduit Flow Process (CFP) for MODFLOW-2005* (Tech. Methods, Book 6, chap. A24). Reston, VA: U.S. Geological Survey.
Spane, F. A., Jr., & Wurstner, S. K. (1993). A computer program for calculating pressure derivatives for use in hydraulic test analysis. *Ground Water, 31*(5), 814-822.

Springer, G. S. (2004). A pipe-based, first approach to modeling closed conduit flow in caves. *Journal of Hydrology, 289*, 178-189. https://doi.org/10.1016/j.jhydrol.2003.11.020

Tsang, C. F., Neretnieks, I., & Tsang, Y. (2015). Hydrologic issues associated with nuclear waste repositories. *Water Resources Research, 51*, 6923-6972. https://doi.org/10.1002/2015WR017641

White, W. B. (2002). Karst hydrology: Recent developments and open questions. *Engineering Geology, 65*, 85-105.

Worthington, S. R. H. (2009). Diagnostic hydrogeologic characteristics of a karst aquifer (Kentucky, USA). *Hydrogeology Journal, 17*, 1665-1678. https://doi.org/10.1007/s10040-009-0489-0

Worthington, S. R. H., & Soley, R. W. N. (2017). Identifying turbulent flow in carbonate aquifers. *Journal of Hydrology, 552*, 70-80. https://doi.org/10.1016/j.jhydrol.2017.06.045