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

Due to the duality in terms of (1) the groundwater ﬂow ﬁeld and (2) the discharge conditions, ﬂow 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 ﬂow, induced by large-scale experiments, in combination with Darcian matrix ﬂow. Effects of large-scale experiments were simulated for idealized settings. Subsequently, diagnostic plots and analyses of different ﬂuxes are applied to interpret differences in the simulated conduit drawdown and general ﬂow patterns. The main focus is set on the question to which extent different conduit ﬂow 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 ﬂow. This work quantiﬁes the relative error that results from assuming laminar conduit ﬂow for the interpretation of a synthetic large-scale pumping test in karst.


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 aly, 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 echal 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 echal 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 5 500 (White, 2002) to N Re 5 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 5 0.1 and f 5 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 acs 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).

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  Exchange flow Q ex [L 3 T 21 ] between the pipe network and the matrix continuum is calculated by a linear quasi steady state exchange coefficient a ex [L 2 T 21 ] as (Barenblatt et al., 1960;Bauer et al., 2003;Shoemaker et al., 2008a) with h c the conduit head [L]. Bauer et al. (2003) introduced the exchange coefficient a ex as a lumped conductance term representing the hydraulic characteristics as well as the geometry of the interface between conduit and matrix as a ex 52pr w Dl p aKs (6) with r w the pipe radius [L] and a a factor which might be interpreted as inverse fissure spacing [L 21 ].

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 with s D the dimensionless drawdown, s the drawdown [L], T the matrix transmissivity [L 2 T 21 ], C D the dimensionless wellbore storage, C the wellbore storage constant [L 2 ], 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 0 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 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.

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.

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 calculated for a conduit length of Dl 5 1 m. A defined signal of Q p 5 0.25 m 3 s 21 close to the pumping well (even distribution of an abstraction rate Q abs 5 0.5 m 3 s 21 in both directions) is used as flow rate according to equations (3) and (4). Assuming a constant flow rate, the head loss differences between laminar and turbulent flow for the conduit diameter of d p 5 2.5 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 5 0.5 m is h c,turb -h c,lam 5 2.5 3 10 23 m whereas a diameter d p 5 2.5 m accounts for head loss differences of h c,turb -h c,lam 5 5.5 3 10 27 m. Even for the same velocity and therefore a higher abstraction rate (Q p 5 6.2 m 3 s 21 ) 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 5 2.5 m. To generate a similar head loss differences of h c,turb -h c,lam 5 2.5 3 10 23 m the flow rate has to be increased to Q p 5 16.7 m 3 s 21 (v 2.5 5 3.4 m s 21 ) with a Reynolds number of Re 5 7.4 3 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 5 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 5 3 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

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 3 113,000 m has a bottom elevation of 0 m with an aquifer thickness of b 5 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 5 1 3 10 24 m s 21 and S s 5 1 3 10 24 m 21 . The well has a constant pumping rate of Q abs 5 0.5 m 3 s 21 , implemented as (negative) Neumann boundary condition (CFPM1 -CRCH Package), for a duration of T total 5 6 3 10 7 s. Groundwater temperature is set to 158C.
All of the following setups consider a single conduit with a length of L p 5 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 aK 5 0.01 s 21 (equation (6)). The roughness height is constant along the conduit with k c 5 1 3 10 24 m, representing a very smooth conduit, and the tortuosity for all conduit segments is set to s 5 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. 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 5 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 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 Dh 1,500 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 3 10 24 m per conduit meter.

Influence on Conduit Flow Pattern: Conduit Diameter
Applying the laminar flow equation on the small diameter conduit (d p 5 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 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 Dh 1,500 m 5 1.8 3 10 22 m. The head loss difference along the conduit between turbulent and laminar flow conditions is Dh c,turb -Dh c,lam 5 9.2 3 10 24 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 5 2.5 m is comparable to that of the laminar flow equation for the small diameter conduit (d p 5 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.

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 5 0.01 m and k c 5 0.25 m (cf., Figure 1b) for a conduit diameter of d p 5 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 5 0.25 m. With a roughness height equal to the radius of the conduit (k c /d p 5 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).

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 3 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 5 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 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.

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 The results presented in Figure 5a can be divided at least into two different parts. Below a relative roughness of k c /d p 5 0.01, the head loss differences are low and uniform. Also the calculated friction factor (based on equation (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 5 0.01, the turbulent friction factor steadily increases. At a relative roughness of k c / d p 5 1, the turbulent friction factor is f 5 0.75 but the head loss difference h c,turb -h c,lam 5 1 3 10 27 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 5 1 is comparable to that for a Reynolds number of Re 5 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 5 0.01 and k c /d p 5 1 describes the transition period between laminar and turbulent flow. For a Reynolds numbers Re 5 10,000 significant head loss differences, able to change the flow pattern, can only be derived by friction factors higher than approximately f 5 10. According to Figure  5b, these friction factors can only be achieved by a relative roughness beyond k c /d p 5 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 5 0.001 m) to the respective conduit diameter (k c /d p 5 1). The conduit diameter ranges between d p 5 0.5 m and d p 5 2.5 m. The values presented in Figures 5c and 5d are based on equations (3) and (4) with a flow rate Q p 5 0.25 m 3 s 21 and no knowledge of the critical Reynolds number and which equation actually applies. The calculated mean velocity ranges between v 0.5 5 1.27 m s 21 and v 2.5 5 0.05 m s 21 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 echal 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 echal 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.

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 5 1 (Jeannin, 2001). According to Figure 5d, those values can only be the result of relative roughness higher than k c /d p 5 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.